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ABSTRACT 


-W@-  consider  the  two-point  boundary  value  problem  for  stiff  systems  of 
ordinary  differential  equations*  For  systems  that  can  be  transformed  to 
essentially  diagonally  dominant  form  with  appropriate  smoothness  conditions,  a 
priori  estimates  are  obtained.  Problems  with  turn^n^  points  can  be  treated 
with  this  theory,  and  we  discuss  this  in  detail.  JWg  give  robust  difference 
approximations  and  present  error  estimates  for  these  schemes.  In  particular 
we  give  a  detailed  description  of  how  to  transform  a  general  system  to 
essentially  diagonally  dominant  form  and  then  stretch  the  independent  variable 
so  that  the  system  will  satisfy  the  correct  smoothness  conditions.  Numerical 
examples  are  presented  for  both  linear  and  nonlinear  problems, 
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NUMERICAL  METHODS  FOR  STIFF  TWO-POINT  BOUNDARY  VALUE  PROBLEMS 


*  **  * 

Heinz-Otto  Kreiss  ,  N.  K.  Nichols  ,  and  David  L.  Brown 

1.  Introduction 

In  this  paper  we  consider  the  two-point  boundary  value  problem  for  a  linear 
system  of  n  ordinary  differential  equations  (ODEs) 

~^-  =  A(x)y  +  F(x) .  OsSx<c,  (1.1) 

OX 

subject  to  n  linearly  independent  boundary  conditions 

B,y(0)  +  Bil/(c)  =g  (1.2) 

Hereyr  =  (y<‘>,  y<2>,  •  ,y<">)  »>  is  a  vector  function  with  n  components  and  B , , 
and  A(x)  e  (?  are  nxn  matrices.  We  assume  furthermore  that  the  vector 
/•(*)  €  C*. 

We  are  particularly  interested  in  the  case  when  the  problem  (l.l)-(l.2)  can 
be  characterized  as  being  stiff  but  not  highly  oscillatory.  The  "stiffness"  of  a 
system  of  ordinary  differential  equations  is  defined  relative  to  a  computational 
grid  on  which  the  system  is  to  be  solved  by  means  of  e.g.  a  difference  approximation : 

ho  Ahi  V-l 

f=^Y  4-B-) — ^ 

x0  *1  *H-1  *N 

We  divide  the  x-axis  into  subintervals  of  variable  length  h.  with  gridpoints 

V 


^ All  numerical  computations  were  made  on  the  Caltech  Applied  Mathematics 
Department  Fluid  Dynamics  VAX. 

Department  of  Applied  Mathematics,  California  Institute  of  Technology  217-50, 
Pasadena,  CA  91125. 

•  • 

Department  of  Mathematics,  University  of  Reading,  England  RG6-2AH. 

•1  If  y  is  a  vector  then  yT  denotes  it*  transpose  and  y*  its  adjoint.  The  vector  norm  is 
defined  by  |y  |  ■  mas  |yW|.  Similar  notations  hold  for  matrices,  for  example 
Ml  *swy  My|/|  ir|.  Furthermore  for  vector  functions  |y(*)|*(  "mas  |y(s)|  denotes 

the  maximum  norm. 

*'  A{x)  c  Cl  if  the  elements  of  A  are  ;  times  continuously  differentiable. 
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Research  partially  supported  by: 

Office  of  Naval  Research  Contract  No.  N00014-83-K-0422;  and  Air  Force  Office 
of  Scientific  Research  Contract  No.  AFOSR-82-0321. 
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S  0,  ty  5  2  ^  ^  =  2| 

j-0 

tions  defined  on  the  grid. 


*n  =  c  and  denote  by  u„  =  u(x„)  vector  func- 
De fining  h  -  max  hj,  we  say  that  the  system  (l.  1)  is 


stiff  if  A  |A|  »  1.  The  case  when  it  is  possible  to  choose  h  such  that  h  |  A  |  «  1 
everywhere  has  been  treated  many  times  before  and  it  is  not  our  aim  to  discuss 
that  situation  here. 

There  are  essentially  two  features  of  stiff  boundary  value  problems  that 
makes  their  solution  by  numerical  methods  difficult.  One  is  that  the  matrix  A 
has  large  eigenvalues  of  both  signs.  The  second  is  that  there  may  be  turning 
points  in  the  problem.  The  concept  of  a  turning  point  is  not  particularly  well 
defined  in  the  literature;  however  for  our  purposes,  we  take  it  to  mean  a  subin¬ 
terval  of  0  <  x  <  c  over  which  an  eigenvalue  of  A  changes  its  order  of  magni¬ 
tude:  In  general,  if  the  eigenvalues  of  A  vary  over  several  orders  of  magnitude 
as  a  function  of  x,  there  are  difficulties. 

It  is  typical  that  the  solutions  of  still  boundary  value  problems  vary  over 
several  different  scales.  Therefore,  it  is  intuitively  clear  that  in  order  to  obtain 
an  accurate  numerical  solution  to  such  problems  the  computational  grid  must 
be  constructed  in  such  a  way  that  the  solution  of  the  problem  is  everywhere 
smooth  with  respect  to  the  grid.  Alternately,  we  can  think  of  changing  the  origi¬ 
nal  problem  by  introducing  a  "stretched”  variable  x  =  x  (x)  such  that  the  solu¬ 
tions  y(2 )  of  the  transformed  problem 


3?'(x)| f-=Ay+F 


varies  slowly  with  respect  to  a  uniform  grid  A„  ■  A.  Since  y(x)  can  change  by 
several  orders  of  magnitude  we  have  to  be  careful  about  what  we  mean  by  "vary¬ 
ing  slowly".  We  must  first  scale  y  correctly.  We  assume  that  this  is  achieved  by 
a  positive  scaling  function  p  e  C*.  The  smoothness  constant  of  this  scaling 


verified-. 


V : 


function  is  defined  by 


K  :=  K\mM\a  :=  sup  max  \(<T<p/  dxv)/  ? 1. 

1  0<Mp  I  • 

/*  A 

Thus  <p  can  grow  like  e*  or  decay  like  s_J&  .  If,  for  example,  A-  =  10,  then  we 
can  obtain  a  change  of  scale  over  a  short  interval. 

Definition  1.1:  A  scaling  function  varies  slowly  for  a  <  z  <  b  with  respect  to 
a  uniform  gridlength  h  if  hK  «  1. 

We  now  consider  y  and  assume  that  we  have  chosen  a  slowly  varying  scaling 
function  <p  such  that  y  -  <py  .  We  can  now  think  of  y  as  being  of  order  0(1)  and 
therefore  we  define  its  smoothness  constant  M'  by 

ft :=  6i„  :=  sup  max  Lf'y  /  dsv  I. 

This  leads  to 

Definition  1.2:  y  varies  slowly  for  a  <  z  &  6  with  respect  to  a  scaling  func¬ 
tion  <p  and  a  uniform  mesh  with  meshsize  h  if 

Kh  «  1,  K.  *  max(K,K) 

The  definition  above  is  useful  for  numerical  purposes  onl>  if  we  can  deter¬ 
mine  the  smoothness  constant  K  without  detailed  knowledge  of  y(x).  In  the 
next  section  we  begin  by  considering  scalar  equations  and  show  that  K  can  be 
determined  in  terms  of  the  properties  of  the  coefficients  of  the  differential  equa¬ 
tion.  In  sections  3  and  4  these  results  are  then  generalized  to  systems  of  ODEs 
which  can  be  smoothly  transformed  to  essentially  diagonally  dominant  form. 

For  the  problem  (1.1), (1.2),  the  constant  K  can  be  determined  in  every 
subinterval  bi*z*btoI0<z*c.  If  K  does  not  change  by  orders  of  magni¬ 
tude  from  subinterval  to  subinterval,  then  we  can  use  a  uniform  meshsize  h  with 
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Kh  «  1  everywhere  on  the  whole  interval.  However,  as  we  pointed  out  above,  it 
is  typical  in  still  problems  that  there  are  several  different  scales  on  which  the 
solution  varies;  for  example  there  can  be  boundary  and  internal  layers.  These 
variations  manifest  themselves  in  values  of  K  which  vary  over  several  orders  of 
magnitude  between  subintervals.  In  sections  9  and  10  we  discuss  how  to  con* 
struct  the  stretched  variable  x  {*)  mentioned  above  so  that  the  new  smoothness 
constant  is  of  the  same  order  everywhere.  This  stretched  variable  leads 

to  a  nonuniform  mesh  in  the  old  variable  x. 

In  the  remainder  of  the  paper  we  discuss  difference  approximations  for  stiff 
systems  and  give  numerical  examples.  There  are  two  basic  classes  of  difference 
approximations  that  can  be  used;  these  are  centered  schemes  and  onesided 
schemes.  Onesided  schemes  have  the  apparent  disadvantage  that  the 
differential  equation  must  be  transformed  to  an  appropriate  form  before  they 
can  be  applied.  However,  centered  schemes  are  not  as  reliable  as  onesided 
schemes,  and  as  we  show  in  section  5,  they  cannot  be  expected  to  give  accept¬ 
able  results  in  general  unless  the  system  has  been  transformed  to  the  proper 
blocked  form.  Even  then,  the  combination  of  onesided  and  centered  schemes 
which  we  advocate  in  that  section  will  perform  better  than  a  strictly  centered 
scheme.  This  is  true  basically  because  the  fundamental  estimates  for  the 
former  more  closely  mimic  those  of  the  differential  equation  than  do  the  esti¬ 
mates  for  the  centered  schemes. 

In  sections  8  through  8  we  discuss  difference  approximations  in  more  detail. 
Section  6  is  concerned  with  difference  approximations  for  scalar  equations,  and 
sections  7  and  B  cover  diagonally  and  essentially  diagonally  dominant  systems. 
In  section  9  we  describe  how  to  transform  a  general  system  to  essentially  diago¬ 
nally  form.  Finally,  in  sections  10  and  11  we  give  a  detailed  description  of  our 
numerical  procedures  and  present  some  numerical  results. 
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2.  Analytic  properties  of  scalar  equations. 

In  this  section  we  consider  scalar  complex  valued  equations 

|j^=  a(x)y  +  f{x), 

y(0)  =  y,.  Oizic  Rea  <  0. 

We  are  interested  in  the  case  that  a(z)  is  large  but  we  ai*e  not  interested  in 
the  case  that  y(x)  is  highly  oscillatory.  Therefore  we  assume  that  there  is  a 
constant  p  ~  0(1)  such  that 

|a/(x)|  < p\aR\  a;  =  Ima,  aR  =  Rea  <  0.  (2.2) 


Later  on  we  want  to  solve  (2.1)  by  difference  approximation  on  a  grid  with 
essentially  uniform  gridsize  h.  Because  our  difference  approximation  will  only 
depend  on  point  values  of  a(z)  and  /(«)  it  is  important  that  the  behavior  of 
these  point  values  represents  well  the  behavior  of  the  continuous  functions.  An 
appropriate  assumption  is  that  there  exist  a  natural  number  p  >  0  and  a  con¬ 
stant  Kx  >  1  with  A, A  «  1  such  that 


l(l«l  ♦1)-‘0lo,*jr..  «<I/H-i)-,|£hI0,**i. 

v  -  1.2... .,p. 


(2.3) 


The  number  p,  which  will  depend  on  the  method  we  use,  will  be  chosen  later. 

We  can  use  (2.3)  to  estimate  the  local  variation  of  a(z)  and  fix).  Letting 

% 

p  -  if  +  sgnf)~ldf  /  dx  and  ^  =  fpsgnf  we  can  write 


2"*  ?(*)/(*)  +  f(*) 


i.e. 


i  «*.r  ll'ution/ 
a/ailaMllty  Codas 
Avail  and/or 
:>ist  i  facial 


»  7  attXf  / 

/(*)*/ •  *  i>irf)dri  ■ft* 
o 


I 


/<*)  =  t"«{,)/<0)  +  *ff8(*) 


where  by  (2.3) 


\g(x)\*Kx,  jy8(x)|  £  K\»  l* 


Therefore  for  |  K\X  |  <  1 

l/(*)-/(0)l  *  l«**‘  - 1|  !/(o)|  +  Ixliail  *  I^1*I(1  +  O(i^i*l)(l/(0)I  +  1). 

A  similar  estimate  holds  for  a(x).  Thus  we  have 
Lemma 2.1:  For  all  x,  with  Kx\xx  -  x0|  <  1 

|o(x,)  -o(x,)|  <  l/f^x,  -x0)|(l  +  0(|#ri(*i  -*0)l))|a(x.)|. 
l/(*i)  -/<*.)!  * \Kx{xx  -x0)|(l  +  0(\Kx(xx  -*0)D)  (|/(x.)|  +  1). 

If  the  expressions  (2.3)  are  not  bounded  then  we  cannot  expect  that  the 
solution  of  (2. 1)  varies  slowly.  As  an  example,  consider  the  differential  equation 


|L=  +  e-K,  y(0)  =  0,  Osxi  1. 


Here  t  >  0  is  a  small  constant  (  for  example  t  =  10  8  ).  Now 


(|a(x)|  +  l)-,|da(x)/dx|  =  |(x  +  2e)"'| 


becomes  arbitrarily  large  for  e  -*  0.  For  xVc  »  1  the  solution  is  to  first  approxi¬ 
mation  given  by 


x  +  e 


and  varies  slowly.  However,  for  0  s  x  *  const.  Ve  it  changes  rapidly.  To  see  this 
we  introduce  a  new  variable  ?  =  x/  Ve  and  obtain 


dy/dx  =  -(z  +  e*)y  +  1,  y(0)  =  0 


Le.  y  (z)  varies  slowly  as  a  function  of  2 '  but  not  as  a  function  of  z.  This  can  also 
be  seen  from  the  graph  of  y(x). 


Pig.  2.1 

The  conditions  (2.2),  (2.3)  still  do  not  guarantee  that  the  solutions  of  the 
differential  equation  change  slowly.  Consider,  for  example, 

dy/dx=-t~ly,  y(0)  =  1,  0  <  e  «  1. 

All  the  above  conditions  are  satisfied.  However  the  solution  is  given  by 

* 

y(z)  =  e“*/* 

which  forms  a  boundary  layer  and  decays  apidly  !,om  1  to  zero.  This  possibility 
can  be  avoided  if  we  assume  also  that  a(0)  is  not  too  large,  i.e.  there  is  a  con¬ 
stant  A*  >  1  with  AgA.  «  1  such  that 


|a*(0)|  <  Km.  is.  )a(0)|  £  (p  +  1)*« 


Observe  that  the  conditions  (2.3),  (2.4)  do  not  prevent  a(x)  from  becoming  large 
because  (2.3)  allows  a  rather  rapid  exponential  growth  of  a(x)  away  from  the 
boundary.  Numerically  this  corresponds  to  the  requirement  that  we  use  an 
exponentially  stretched  mesh  to  represent  boundary  layers  in  the  solution. 

The  above  conditions  are  not  sufficient  to  guarantee  that  the  solution  y(x) 
stays  bounded.  Therefore  we  assume  explicitly  that  the  equation  has  been 


properly  scaled  such  that 


flv(*)!lo.e  <  !• 


The  following  lemma  says  that  (2.5)  implies  the  condition 
1/  I /(I®  I  +  1)  <  constant. 

Lemma  2.2:  Assume  that  c  >  2/  Kx.  Then  there  is  a  constant  C3  which 
depends  only  on  Kx  such  that 


II  |o(x)|  +  1  H°'8  A  C>SIIV  H°,e 


(2.5a) 


Proof.  Without  restriction  we  can  assume  that  Dl/llo.c  =  1-  bet  x0  be  a 


puiut  with 


/(*.) 


If  |  c  -  x,  1  ;e  1/  Kx,  we  consider  (2. 1)  for  x  >  x,  and  introduce  new  variables 

(»  -*»)  ~  _  !»(*•)!  +  1 

I«(x.)|  4  r  V 
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dy/dx  =(1  +  |a(x,)|)-‘B(*)y  +/(*)/ |/(*,)| . 


By  lemma  2.1 


1  0  (x  )  I 

a(x„) 

+ 

a(z)  -  a (x0 ) 

1 a  (x, )  |  +  1 

|a(x,)|  +  1 

|a(x0)|  +  1 

*  1  +  0(1  *,(* -*.)!). 


and 


l/GOI 


^  1  —  0(1  Ki(x  -  *0 )  I  )■  Therefore 


implies 


II V  ll*#.,,i  I V  (*  i)  I  *  (*i  ~*o)(l  -  0{\Ki(xi  -*0)l)  - 

(x,  -*0)(1  +  0{\Kl(zl  -z0)|)||gr  ||,#^. 


Thus  for  |x,  —  x0  |  =  6/  Kx.  S  >  0  sufficiently  small. 


|o(*o)l  +  1 

!/(*.)! 


\v  I 


*8-»l' 


Since  the  last  inequality  holds  for  any  z, ,  the  inequality  of  the  lemma  holds  in 
this  case  with  C3  =  2X\/<5.  If  |c  -  x,  |  <  1/ Kx  then  by  considering  (2.1)  for 
x  <  x,  instead,  we  obtain  the  desired  inequality  in  a  similar  way.  This  proves  the 
lemma. 

We  shall  now  prove  that  under  the  conditions  of  (2.2)-(2.5)  there  is  a  con¬ 
stant  C  such  that  the  solution  y(x)  of  (2.1)  satisfies 


L 


||dwy/olx,'||o,e  s  C  . 


v  =  1,2,..., p. 


(2.6) 


This  is  easiest  to  show  when  |a(z)|  is  small  tor  all  z:  We  have 


Lemma  2.3:  Assume  that  (2.5a)  holds,  that  conditions  (2.2)-(2.5)  are 
satisfied  and  that 

l|a*(*)llo.e  * 

Then  there  is  a  constant  C  which  depends  only  on  AT,,  Af2,p  such  that  (2.6)  holds. 
Proof:  (2.2)  implies 

II  “(*)  llo.e  ^  (P  +  1)*2. 

Therefore  by  (2.3) 

|| dva/ dxv\ 0-c  s  /r,((p  +  l)Ag4-l)  i/  =  1.2 . p. 

By  (2.5a)  and  (2.3) 

ll/(*)llo.e  *  Cs(l|tt|lo.e  +  1)  ||  dv/ /<**’' II  o.c  <  const..  1/ =  l,2....,j>. 

Using  the  differential  equation,  (2.6)  follows. 

We  assume  now  that 

oB(z)  <  -l.  (2.7) 

i.e.  we  allow  |  a#  |  to  become  arbitrarily  large.  We  will  use 

Lemma  2.4:  The  solution  of  (2.1)  with  a#(z)  =  Rea  <  0  satisfies  the  esti¬ 
mates 


|y(z)|  £  z  max  | /“ (17)}  +  s 

0<T)<X 


7  •«(««< 


ly  (0)  I 


(2.8a) 


|y(z)|  *  max  1/ (tj)/ aJ?(z)|  +  e 

0«y<s 


m 

f 


lv(o)l 


(2.8b) 
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Proof:  The  solution  of  (2.1)  can  be  written  down  explicitly: 


a  a 

*  /•(<)*<  /•(«*( 

v(*)  -  /  e”  f(v)^v  +  «°  v(0). 

o 


The  first  estimate  follows  from  the  inequality 


a  a 

f  •«)<*«  /  •*(()<*( 

|e  *  |  «  e  ’  i  1, 


Furthermore 


*  7  «*«)•<<  7  «*«)<< f 

\y{x)\^fe’>  \f{ri)\dr)  +  |e  0  I  ly(0)|  < 

o 


*  /«*(<>*< 

7  aR^rj)e  v  \/(v)/*Jt(v)\*V 
o 


f  «*«)<** 

+  e°  Jy(0)|  -C 


d  /•»(«-« 

J  7^re  dr) 


/■*(«<« 

max  |/(t7)/a*(»?)l  +  «°  ly(0)l 

0*r^a 


which  gives  us  (2.8b). 

Using  lemma  2.4  we  can  obtain 

Lemma  2.5:  Consider  (2.1)  and  assume  that  the  condition  (2.2)  is  satisfied. 
Then  there  is  a  constant  Cx  which  depends  only  on  p  such  that 

\dvy(x)/ dxv\  s  (2.9a) 


Cx 


v-\ 

£ 

f»0 


♦i. 


lo* 


^(0) 

dxv 


,v  =  0,1.2 . p. 


If  furthermore  (2.3),  (2.5a),  and  (2.7)  hold,  then  there  is  a  constant  Cj  which 
depends  only  on  Kx  and  p  such  that 


* 

•i 


i 
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\dvy(x)/dxv\  <  Ct  Uy  ||  o.»  + 


(2.9b) 


Proof:  For  1/  =  0  the  first  two  estimates  follow  from  (2.2)  fend  lemma  2.4 
since 

!//«**!  s  l//«*lla/a#l  *  (1  +  />)  l//«l 

Now  let  u  =  dy/dx.  Differentiating  (2.1)  gives  us 

§■=““  +  (f^  *  (210> 

Thus  by  lemma  2.4 

\\du/dx\\0'M  <  II  a/? 'da  /  dx\\  0-J,  |j  y  ||  <>.#  +  l|a*  ld/ /  dx  ||0-,  +  |dy/dx|,M0. 

Now 


\\aj'da/dx  ||0iJ,  i  (1  +  p) || o'^o/ d*  ||0i, 

||  a* 'd/  /  dx  ||  o.,  *  ( 1  +  p)  ||  a"'d/  /  dx  ||  0jt  s  ( 1  +  p)  ||  II  II  o., 

and  if  |  a  |  Js  1  we  have  also  that 


and 


||o-,da/  dz  ||q^  <  2|j(|  a  |  »  l)~’do/dx  ||0i. 


I  (1/  L»  g  2|  »  ! 


and  so  (2.9a,b)  follow  for  u  =  1  from  (2.2),  (2.3),  (2.5a)  and  (2.7). 

The  estimates  for  higher  derivatives  are  obtained  by  repeated 
differentiation  of  (2. 1).  This  proves  the  lemma. 


'  « 
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We  can  now  prove  the  main  result  of  this  section 

Theorem  2.1:  Assume  that  C*2/Kt  and  that  the  conditions  (2.2)-(2.5) 
hold.  Then  there  is  a  constant  C  which  depends  only  on  Kv  Kt.  p  each  that  (2.6) 
is  valid. 

Proof:  Lemma  2.2  tells  us  that  (2.5a)  holds.  Now  divide  the  interval 
Oiz^c  into  as  few  subintervals  cr<z  £  cril  as  possible  such  that  at  least  one 
of  the  two  conditions 

holds  in  the  whole  subinterval.  If  J  a* )  £  Kg  is  valid  then  we  can  estimate  the 
derivatives  by  lemma  2.3.  If  a#  <  —  1  then  we  can  obtain  the  estimate  using  lem¬ 
mas  2.2  and  2.5  if  we  have  a  bound  for  £  \d*y(cr)/ dx*  | .  The  interval 

i-o 

0  =  c0<zs;c1is  included  in  the  case  |  a*  |  «  Kg.  so  we  are  only  concerned  with 

the  remaining  intermediate  points  t  =  1,2 . For  these  points  we  obtain  this 

bound  from  the  estimate  for  the  previous  subinterval  cT_,  <  x  <  cr  This  proves 
the  theorem. 

Finally,  we  can  eliminate  the  condition  (2.5),  i.e.  the  assumption  that  we 
have  scaled  the  solution  beforehand.  If  fly  ||q>6  >  1  then  y  =  y/  Hyllo.t  satisfies 

dy/dx=a(x)y  +  f ,  f  =//!vlio.e  (210) 

Now 

8(17  I  \)~'dvf  /dz*'|lo.e  s  11(1/  I  +  l)-‘dV/d*i lo.. 
implies  that  (2.10)  satisfies  all  our  conditions.  Thus 
Hd^/dz-Ho.,  *  C, 

i.e.  ||  d1'!/ /  d**' ||  0<e  s  C || y  ||  e,e  for  H  V  llo.c  >  1- 


(211) 


Combining  (2.6)  and  (2.11)  gives  us 

Theorem  2.2:  If  the  conditions  (2.2)-(2.4)  are  satisfied,  then  we  obtain  the 
estimate 


S«f*'v/d*‘'llo.e  <  C(||y  it 0.e  +1).  V-  1.2 . p. 


(2.12) 


Remark:  If  hC  «  1  then  y(a)  is  slowly  varying  with* respect  to  the  scaling 
function  9  =  9  y  ||  o.e  •  We  can  obtain  a  much  more  precise  result  by  using 
p  =  V(|/  I*  +  1)/(| a  |*  +  1)  as  a  scaling  function.  To  explain  this  we  assume 
that  a  »  1,  /  »  1  and  use  9  -  |  /  /  a  (.  Then 


dp/ da 
<P 


*  l/V/l  +  | a’/ a |  ~2tf 


and 9  =  y/p  satisfies 


fjf*  s  (a  -  9'/9)V  +  “•  <•••  I# I  “  |a/(a-p)|  ~  1. 
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3.  Diagonally  dominant  systems. 

We  consider  now  systems  of  differential  equations.  A  reasonable  assumption 
for  such  systems  is  that  the  coefficients  change  slowly,  Le.  that  they  satisfy  con* 
ditions  of  the  same  type  as  we  described  for  scalar  equations.  Unfortunately, 
this  assumption  is  not  sufficient  to  guarantee  that  the  solutions  of  the  systems 
also  vary  slowly.  Consider,  for  example,  the  system 


d_* 
dx  v 


i/(-l)  =  a.  v(l)=/?  (3.1) 

In  a  neighborhood  of  x  =  0  the  functions  x*.  e"1,  e*x  satisfy  the  conditions  (2.3) 
with  Kx  =  1.  However,  the  solution  is  not  smooth.  In  component  form  (3.1)  is 
given  by 

u'  =  xhj ,  ev’  =  u  -  c*x 
i.e. 


ev" 


=  z*v  -  e*. 


A  graph  of  the  solution  is  given  in  fig.  3.1. 


There  is  however  a  class  of  problems  that  behave  like  scalar  equations  and 
we  shall  describe  this  class  now. 

» 

Definition  3.1:  The  matrix  function 


A(x)  = 


“!•»(*) 

oa,(x)  • 

Vi(*)  •  •  •  • 

"  „  *  V 


is  diagonally  dominated  if 
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Rea*  <0  <  =  1.2,. ...r. 

Rea«  >0  <  =  r  +  1 m. 

ft)  There  is  a  constant  6  >  0  independent  of  z  such  that 

§  I  36  (l  —  <5)1  Rea*  J,  i  = 
f» H 

7)  There  is  a  constant  p  ~  0(1)  such  that 

|lma«|  <  p|  Rea«  | .  i  =  1,2 . m, 

If  A  is  diagonally  dominated  then  it  is  appropriate  to  write  the  system  (1.1) 
in  the  form 


^-  =  A(/  +  5)y  +  F,  Oiiic 


(3.2) 


where 


fi1  =  (*00(0,1.022.  •  •  •  lOft).  fia  =  *ag(aTM.rM.<lr**.r*r  •  •  • 


Such  a  system  is  said  to  be  in  diagonally  dominant  form.  Corresponding  to  (2.3) 
and  (2.4)  we  make 

Assumption  3.1:  There  is  a  constant  Kx  such  that 


I  A"'  —■#  o,  *  AT,.  ||  ~f-||  0.e  <KX.  v  =  1.2,..., p  . 


Assumption  3.2 


l|A-1“~^O.eS^1.  «/  =  1.2..... p. 

ax 


Romark:  Corresponding  to  the  previous  section  we  can  replace  assumption 


3.2  with 


We  shall  always  assume  that  |a«|  »  1  and  therefore  assumption  3.2  means  that 
we  have  already  scaled  the  equation  properly. 

Assumption  3.3:  There  is  a  constant  Ka  such  that 
^(0)1**,.  |A"(c)|s  AT,. 


We  shall  show  that  under  these  conditions  the  solutions  of  the  system  (3.2) 
change  slowly.  We  start  again  with  a  couple  of  lemmata: 

Lomma  3.1:  Assume  that  A(z)  is  diagonally  dominated.  Then  the  solutions 
of  (3. 1)  satisfy  the  estimate 


(3.3) 


lv(*)l  <  l4l!(A  +  h')~lFh.»  +  *(*)] 

M  \  4 


where 


Hdj  •/U)*t  M»  J  •"«)*€ 

*(*)  =  f  9  lv'(*)l+*  *  lvff(c)|.  (3.4) 


Here  we  have  used  the  notation 


y*(x)  =  (y{,).  •  •  •  .ytr))r.  ytf(*)  =  (ytrM).  •  •  •  .ytm))r 

o/(x)  =  min|ReoM(*)|,  a//(x)  =  min  |Rea*(x)|  . 


Proof:  We  consider  first  the  case  y/(0)  *  y;/(c )  =  0.  Let  M  denote  the  space 
of  all  continuous  functons  g  with  gl( 0)  =  gn(c)  =  0.  The  differential  equation 


Uv  »  |j[--Ay  */* 


has  a  unique  solution  in  ■  given  by 


=  U'F  = 


•  ftfmt 
/• 9  F‘(rl)  dr, 

0 

[/>  ^(q)di7 


Let  L|  denote  the  operator  defined  by 

L,V  =  {A  -  A)y,  y  €  M. 

Then  we  can  write  the  differential  equation  (3.2)  in  the  form 

(/  -  V%)V  =  W'F 


Is  the  same  way  as  in  lemma  2.2  it  follows  that 
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HVl/(*)Ho.e  *  2||(A+  AV/’OOIIo.. 


Furthermore 


I  V'LtV  II o.e  *2|](A  +  A*)-‘L,V  Ho..  ^  (1  "  «>llv  Ho.« 


and  the  estimate  (3.3)  follows  for  the  case  that  y7(0)  =  yn(c)  -  0. 

Now  we  consider  the  case  where  /  *  0  and  yl(0),  yri(c)  are  arbitrary.  We 
let  y  be  the  solution  of  (3.2)  and  write 


_  fy'(o)1  • 

y  =  v  +  Dy„  y,  =  D  - 


m 


0  e 


Then  v  satisfies 


V  =  Au  +  (A  -  A )Dy„  v7(0)  =  v/7(c )  =  0. 


By  (3.3)  we  obtain 


M*)io.e  *  i_,»2(A  +  A')"*U  -  A) II o,.  Iloilo..  •  Iv.  I  * 

* 4"  ■  n,“l?,i^rj lv>  1  *  s"(l  • 61  '»• 1  - 


I v (* ) I  *  -  i)  ly.  I  +  Iv.  I  =  i"‘ly.  I- 


To  show  that  in  addition  to  being  bounded,  y(x)  decreases  exponentially  away 
from  the  boundary,  we  consider  the  solution  of  (3.2)  with 


/  ■  0,  y7(0)  =  y7,  y77(c )  =  0. 


Let  y  =  ,  p(x)  =  o 7(()<*(-  Then  *(x)  satisfies  the  equation 


-20- 


«•(*)  =  U(x)  «•  f»'(* )/)*(*). 

Using  (3.5)  we  obtain 

Therefore 

|y(z)|  £e*<*>|«(*)l  <8  2«'1e*W|y/l- 

The  corresponding  result  holds  for  yr(0)  =  0,  yw(c)  =  y[l.  Thus  we  have  proved 
the  complete  estimate  (3.3). 

In  the  same  way  as  for  the  scalar  equation  we  can  now  use  the  last  lemma 
to  discuss  the  smoothness  of  the  solutions  of  (3. 2). 

Lemma  3.2:  Consider  the  system  (3.2).  Assume  that  A(x)  is  diagonally 
dominated  and  that  the  conditions  of  assumption  3.1  are  satisfied.  Then  there  is 
a  constant  Ct  which  depends  only  on  AT,,  p,  and  6  such  that 


Ptoof:  The  proof  resembles  the  pioof  of  lemma  2.5  elosuly.  Foi  u  =  0  the 
estimate  is  given  by  lemma  3.1.  Let  u  =  dy/  dx  and  differentiate  (3.2).  Then 


=  A(z)u  +  F.  + 

dx  dx 


dA 

dx 


(3.8) 


Using  the  estimate  for  y  we  obtain  an  estimate  for  F  '  and  lemma  3.1  gives  us 
the  estimate  for  u .  This  process  can  be  continued  and  the  lemma  is  proved. 

We  can  now  prove  the  main  result  of  this  section. 


-21- 


Theorem  3.1:  Consider  the  diagonally  dominant  system  (3.2)  and  assume 
that  assumptions  3. 1-3.3  hold.  Then  there  is  a  constant  C  which  depends  only  on 
Kx,  Ri,  p,  and  S  such  that 

8  ^Ho,  <  +  l].  v  =  1.2... ..p.  (3.7) 

Proof:  By  lemma  3.2  and  assumption  3.2 


By  assumptions  3.2  and  3.3  and  since  A  is  diagonally  dominated, 

£  lay  |  <  2Kt.  |/<*>|  <  K2RX.  t  =  1.2 . r 

J-i 

holds  at  *  =0.  Therefore  the  differential  equations  gives  us 

<  2A*|y(0)|  +  K2RX. 

Correspondingly  we  get  for  x  =  c 

*2^  s  2/fg|y(c)|  t  K2RX. 

Thus  by  lemma  3.2  we  can  estimate  (( dy/  dx  ||0^.  The  estimates  for  higer  deriva¬ 
tives  are  obtained  as  before  by  differentiation. 

If  A7(0),  A rr(c)  are  not  0(1)  then  we  introduce  an  exponential  stretching 
such  that  they  are  bounded  in  the  stretched  variables.  Let 

a,  =  max  |  a*  (0)|  >  1 


and  introduce  new  variables  by 
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x  =  of1/ *(()<*$•  *(()  =  •  0 
0 

where  g  will  be  defined  below.  Then  the  system  (1.1)  becomes 


Jjf-=  a  f  *  p(z  )Ay  +  a“V(2)/,(*).  (3-8) 

Let  x  i  =  /f,-,loga,  and  choose 


g(v)  = 


AT,  for  0  <  t)  <  x  ! 
0  for  r)>x  i 


is. 


*(2)  = 


K.tt  . 

i  1  for 
(a,  for  x 


Oil <z , 

>*i 


Then 


z 


-  lJ/aj/Tj  for  Oil  Si, 
?+*!-*!  for  I  >*! 


z,  =  (1  -o  f‘)/AV 


Now  treat  the  neighborhood  of  z  =  c  correspondingly.  Then  A/(0),  A/7(c)  are 
bounded  and  assumptions  3.1  and  3.2  hold  for  p  =  1  with  Klt  replaced  by  2 Kx, 

2 /tv  Thus  the  estimate  of  theorem  3.1  is  valid  for  p  =  1.  In  particular  || 

(where  xt  =  c  |1  aa*’ l/ATi,  and  ag  =  max|a*|)  is  already  bounded  in  the 

unstretched  variable  because  for  z,  <  x  <  zK  no  stretching  occured. 

To  obtain  estimates  for  higher  derivatives  we  could  replace  p(rj)  by  a 
smoother  function.  However,  this  is  not  necessary.  Apply  the  stretching  to  the 
differential  equation  (3.6)  for  u  =  dy/dx.  The  we  get  a  bound  for  du/dS.  On 
any  of  the  subintervals  0  £  z  £  zt,  zt  £  z  £  Zg,  zg  £  z  £  c  differentiation  com¬ 
mutes  with  stretching  and  therefore  we  can  estimate  the  derivatives  on  every 
subinterval  in  terms  of  j|y!lo.c  +  1.  In  particular  we  have  for  the  subinterval 
away  from  the  boundaries: 


Theorem  3.2:  Assume  that  assumptions  3.1  and  3.2  hold  and  that  c  %  2/  Kx. 
Then  there  is  a  constant  C  which  depends  only  on  K\,  Ag.  p  and  6  such  that 
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4.  Essentially  diagonally  dominant  systems. 

The  class  of  diagonally  dominant  systems  is  not  broad  enough  to  include 
many  interesting  problems.  Therefore,  in  this  section  we  generalize  our  results 
to  systems  which  are  essentially  diagonally  dominant  or  which  can  be 
transformed  smoothly  to  systems  of  that  type. 

Definition  4.1  A  matrix 


“11 

0,2 

°ai 

022  ’ 

aa» 

°ni 

•  - 

is  called  essentially  diagonally  dominated  if  there  is  a  constant  Ko  with  K„h  «  \ 
such  that  4  can  be  partitioned  into  the  form 


A 


4,2 


Au  =  A (/  +  B) . 


0 

ass 


where  A,,  is  diagonally  dominated  and 

II A  *4,2 1|  0e  ^  K  •  IMyllox  SA,,  j  =  1,2. 


We  consider  now  systems  ( 1.1) 


(4.1) 


where  4  is  essentially  diagonally  dominated.  Such  a  system  is  said  to  be  in 
essentially  diagonally  dominant  form.  We  are  interested  in  the  case  that 
A  |  A,,  |  »  1.  We  make 


K 


Assumption  4.1: 
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HA-'d^A/dill,..  *  AT,.  || bTxdvAxi/  d*l'll0.e  «  Kt. 

UvB/dxv\ |0.e  *  Kt.  \\dvAZi/dxv || o.e  *  Kv  j  =  1.2;  i/  =  1.2,....p  . 

Assumption  4.2: 

||A-‘d‘'C/iir‘/||o.c  ^  ^i.  \\dvH/dxv\\0ft  zXt.  v  =  0.1,2... .p. 

Assumption  4.3: 

|A'(0)|*tfa.  \\a(c)\*Kt. 

We  want  to  show  that  the  estimate  of  theorem  3.1  is  still  valid. 

Theorem  4.1:  There  is  a  constant  (?  which  depends  only  on  Klt  Kz,  p  and 
6  such  that 

\\dvy  /  dxv\\fix  <  ?<||y  ||0.«  +  l).  V  =  (4  2) 

Proof:  We  have  by  (4.  l)  and  assumptions  (4.1)  and  (4.2)  that 

||  dtu  /  dx  1  o,c  ss  ATjy  Ho.e  +  %i- 
Write  the  equations  for  y  in  the  form 

dy/dx  =  Auy  +  +  G 

where  Gj  =  Alzw .  Now, 

|!  A-'cfGj/  dx  ||  0,e  ^  K.  ||  dw  /  dx  ||  0>c  +  Kjw  ||  0,e  • 

so  in  the  same  way  as  the  proof  of  theorem  3.1  we  obtain 


\\dy/  dx  ||  0  e  £  const.  (||y  ||e>t  +  1)  ■ 
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This  proves  the  theorem  for  p  =  1.  For  higher  derivatives  it  follows  by 
differentiation  of  the  differential  equations. 

If  the  assumption  4.3  is  not  satisfied  then  we  obtain  from  theorem  3.2  that 
the  derivatives  are  still  bounded  in  the  interior  of  0  s  x  £  c.  Thus  we  have 

Theorem  4.2:  Assume  that  in  the  neighborhood  |*,-*|  <  p,  pit  2/  Kx  of 
every  point  x9  with  1/  K\  <  x9  sfi  c  -  1  /  Kx  we  can  write 'the  system  (1.1)  in  the 
form  (4.1)  such  that  assumptions  4.1  and  4.2  hold.  Then  we  obtain 

H*VV  /  <**1  i/Kl.e-l/Kl  S  C(I  V  II  0..  +  1).  (4.3) 

Remarks:  It  is  important  to  note  that  the  dimension,  m,,  of  the  large  block 
need  not  be  a  constant  on  the  entire  interval  0  <  *  <  c.  As  the  assumptions  of 
the  theorem  state,  it  is  only  necessary  that  we  be  able  to  block  the  system  into 
the  form  (4.1)  in  the  neighborhood  of  every  point  in  the  interior  of  the  interval. 
Thus  the  theorem  applies  to  problems  with  "turning  points",  since  by  this  we 
mean  a  problem  in  which  one  of  the  eigenvalues  of  the  large  block  <4lt  locally 
changes  size  by  an  order  of  magnitude. 

Observe  also  that  the  estimate  (4.3)  is  invariant  with  respect  to  smooth 
transformations,  i.e.  we  can  replace  y  locally  by  y  =  S(x)y  where 

151*-***+  US-1  II 

Uvs/dx'’}\, Kx,  v  -  1,2, ...p,  <4,4) 

and  an  estimate  of  the  type  (4.3)  is  still  valid.  Therefore  we  can  relax  the  condi¬ 
tions  of  theorem  4.2.  Instead  of  assuming  that  A  has  the  form  (4.1)  we  need  only 
assume  that  in  the  neighborhood  of  every  point  there  is  a  transformation  S . 
satisfying  (4.4).  such  that  SAS~l  is  of  the  form  (4.1)  and  satisfies  the  conditions 


of  theorem  4.2. 
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We  shall  now  derive  estimates  for  essentially  diagonally  dominant  systems. 
We  consider  An  to  be  the  large  part  and  Aig,  Ag |.  Au  the  0(1)  part  of  the  system 

(4  1). 

It  is  well  known  that  the  solutions  of  the  system 

~~  -  AggW  +  F,  0  <  x  «  c 
02 

satisfy  an  estimate 

Mo*  ^(MO)!  +c||^||0.e).  Afa  *ejep[||Agg||0.e  <?]•  (4.5) 

Note  that  since  we  do  not  assume  that  the  eigenvalues  k  of  Agg  satisfy  Re«c  <  0, 
this  bound  for  Kg  is  realistic  and  so  the  estimate  (4.5)  is  useful  only  if  c  is 
sufficiently  small. 

We  can  now  estimate  the  solutions  of  (4.1)  in  terms  of  G.  H  and  y7(0), 
y7/(0),  to(0).  Lemma  3.1  and  (4.5)  give  us 

Bvilo.c^2d  1 II  (A  +  A  )_,i4jg||ofe  H'"  II  o.c  +  D\  • 

8 w  II  o.e  <  K gc  ||  Ag\  ||  o.c  II  v  II  o.c  +  Dg  , 

where  Dx  =  2d-‘(l|(A  +  A*)“lG||0>e  +  |y7(0)|  +  |y7/(c)|).  and 

Dg  -  /Ti(|u»(0)  |  +  c  fl//Ho.c).  Therefore  we  have 
Lemma  4.1  Assume  that 

2d~‘cA'88(A+  A*)'li4i8lo^  Mg,ll<u*  1-d*.  6*  >  0. 

then 

*  *  I V  8  o..  *  Si  +  2 d-»  ||  (A  +  A*)_,Aig||  o.«  Dg, 

<**Mo.«  4  Dg  +  Kgc  ||Agillo.e0i. 

i 


(4.8) 
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if  c  is  not  small  then  we  can  derive  global  estimates  in  the  following  way. 

We  divide  the  interval  0  <  x  <  c  into  subintervals  c4  <  *  s e  c4M,  i  =  0.1.2 . g-1, 

c,  =  0.  cf  =c.  c4f,-c4  sufficiently  small.  On  every  subinterval  we  write 
y  =  yF  +.  yH,  w  =  t up  +  w/f  where  y P,  % up  denote  the  solution  of  the  differential 
equation  with  boundary  conditions 

yfai)  =  V^(c4*i)  =  tuj>(c4)  =  0.  (4.7) 

and  yff,  to#  is  the  solution  of  the  homogeneous  differential  equation  with 

yk(ci)  =  yf(c4).  yflicm)  =  y//(c4«.,).  wH(c{)  =  so(c4).  (4.8) 

By  lemma  4.1  yf>(c 4M).  and  u»/»(c4*j)  are  bounded  and  are  uniquely  deter* 

mined  by  (4.7).  Also,  remembering  that  the  differential  equations  ore  linear,  we 
can  write 


vff(Ci)  =  P".  yk(Ci* i)  =  «"*(*<♦ i>  =  P 

Here  Pn ,  P1.  and  P  are  linear  relations  with  bounded  coefficients  in  y/(ct), 
y7/(c4+l)  and  tu(c4).  Thus  in  every  subinterval  we  obtain  n  linear  relations 

yir(ci)  =  Pa  +  yff(ci).  yl(cul)  =  Pl  +  y£(c4<.|).  tu(c4M)  =  P  +  -u>j»(ct M) 

for  the  variables  y (c4).  y(c44.i).  ui(c4),  and  lotc^).  There  are  n(g  +1)  unknowns 
y(cj),  tu(cj)  and  g  subintervals.  The  missing  relations  are  obtained  from  the 
boundary  conditions  for  the  original  problem.  Thus  the  y(cj),  w(cj)  can  be 
obtained  as  the  solution  of  a  linear  system  of  equations.  Whether  we  obtain  rea* 
sonable  bounds  for  the  original  problem  (1.1).  (1.2)  depends  on  the  condition 
number  of  that  linear  system. 


.  V  *  -  .y  /K'.; 


4 


1 
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6.  Tine  choice  of  difference  methods 

In  this  section  we  shall  discuss  our  choice  of  difference  methods.  Perhaps 
the  simplest  stiff  problem  is  given  by 

-V  ■*•/(*)  •  V(°)  =  V*  (5.1) 

where  0  <  e  «  1  is  a  very  small  positive  constant  and  /  (* )  is  a  smooth  function 
with  derivatives  which  are  0(  1).  The  solution  y  s  ys  +  VB  consists  of  a  smooth 
part 


y s(*)  s  /(*)  +  0(e). 


(5.2) 


which  can  be  obtained  by  an  asymptotic  expansion,  and  a  boundary  layer  part 


VB  -  «  (v.  '  Vs(0)). 


(5.3) 


which  is  a  solution  of  the  homogeneous  equation 


(5.4) 


Thus  the  solution  is  smooth  except  in  a  boundary  layer  near  x  =  0  where  it 
changes  rapidly. 

Now  consider  a  uniform  grid  x„  =  vh,  v  =  0,1,2,...;  0  <  h  «  1.  There  are  two 
standard  types  of  difference  approximations.  One  type  is  centered  schemes  of 
which  the  trapezoidal  rule  is  an  example: 


4  U„  _  /  y+l  + 

* - - - - - - - 4- - r - .  u,  -  y. 


(5.5) 


The  other  type  is  one-sided  schemes,  such  as  the  implicit  Euler  method: 


v„M  - 

t - r - s  +  fv*\'  V.  ~V% 


(5.5) 
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The  solutions  u„  =  ti$„  ♦  ug^  vv  =  Vgv  ■*■  Vgv  of  (5.5)  and  (5.6)  respectively,  con¬ 
sist  again  of  a  smooth  part 

=  /*+  0(e).  v¥  a  vSv  +  0(e) 

and  a  boundary  layer  part 

«*  \  ■  <5-7) 
«i.  =  T"(v.  -VS.).  r  =  ytV/T 

where  /cv.  t*'  are  solutions  of  the  corresponding  homogeneous  difference  equa¬ 
tions.  If  h  «.  e  then  k  ~  t  ~  a~h/t  and  therefore  kv  ~  r*  ~  i.e.  the  solutions 

of  the  difference  schemes  approximate  the  solution  of  the  differential  equation 
well.  However,  we  are  interested  in  the  case  that  e  «  h.  In  this  case 

*~-l.  |t!~0.  (5.7a) 

Thus  u¥  is  in  general  highly  oscillatory  everywhere  and  does  not  approximate 
y(x)  well  at  all.  In  contrast,  \vv  -  y(xv)  |  is  small  away  from  the  boundary  layer. 
The  advantage  of  one-sided  methods  in  this  situation  is  clear. 

Onesided  schemes  have  a  major  drawback,  however,  when  they  are  to  be 
used  for  solving  systems  of  equations.  For  the  equation 

Cdt~V  *  f '  V(c)=V..  x*c, 

the  appropriate  onesided  scheme  is  the  explicit  Euler  method, 

tV~~h  Vv  ~  v\>  +  /»■  *•••  +  vjv  =  Vo  .  (5-8) 

because  we  start  the  integration  at  z  =  c  and  calculate  the  solution  for 
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decreasing  values  of  v.  This  construction  of  onesided  schemes  can  be  general¬ 
ized  to  systems  of  the  form 


dy 

dx 


ki  o  o 

V1 

° 

^  ° 

o  © 

y  +  By  +  F, 

V  = 

yn 

VIU 

(5.9) 


where 


|/L4o|  «  1.  \hB |  «  1.  | |  »  1. 


and  the  eigenvalues  k(A- i),  cCd+i)  of  A-i  .  A+\  respectively  satisfy  the  inequali 
ties 


-h Re  *c(A_t)  »  1.  hRe  jc(A+i)  »  1. 

An  approximation  to  (5.9)  on  a  nonuniform  mesh  is 

*oi„. 

fly 

u".'-u".  =  ^  (2v)u//  +  ^  (Xwvl)u//v|)  +  +  c"  ! ).  (5. 10) 

,,///  _  ..III 

Aw 


where 


G„  =  B{x  w)uK  +  /’w. 

t.e.  we  use  implicit  or  explicit  Euler  for  the  variables  corresponding  to  "large” 
eigenvalues  of  negative  or  positive  sign,  respectively,  and  we  use  the  trapezoidal 
rule  for  the  variables  corresponding  to  "small"  eigenvalues.  If  the  system  (l.l) 
is  not  already  of  the  form  (5.9),  however,  we  must  transform  it  to  that  form 
before  we  can  tell  which  combination  of  these  methods  to  use.  Since  this 


i  -*  wmrf 

SiWifltfiMii-i  n  rfitt't  dVuliil  in  «*i 
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transformation  can  be  expected  to  be  somewhat  difficult  to  implement  numeri¬ 
cally,  the  question  arises  as  to  whether  centered  schemes  can  still  be  success¬ 
fully  employed  for  stiff  problems,  since  they  do  not  require  a  priori  knowledge  of 
such  a  transformation  before  they  can  be  written  down.  The  answer  is  that  in 
many  cases,  centered  schemes  can  be  used  together  with  appropriately  chosen 
nonuniform  m§sh*s.  Consider  again  the  trapezoidal  rule  (5.5).  Instead  of  using 
a  uniform  mesh  we  now  use  a  nonuniform  mesh  made  up  of  two  uniform  meshes 
with  meshwidths  h  «  e  and  h  respectively. 


h  h 


In  the  boundary  layer  0<x  s*.  X  -  O(e|loge|)  we  use  h  and  return  to  h  for 
x  >  X.  In  this  case  the  boundary  layer  part  of  the  solution  is  given  by 


~«So) 

±(-l  -Mfc)  Xv>£ 


It  is  clear  that  by  choosing  £  sufficiently  large  we  can  make  \uSv  -  ys(x„)  I 
small  everywhere.  For  systems  we  proceed  correspondingly.  We  use  a  fine  grid 

i 

in  the  neighborhood  of  x  =  O.c  and  a  coarse  grid  in  the  interior.  On  this  coupled 
grid  we  approximate  (1. 1)  by 


-4"""=  fc(i4  (*„♦!)«.»♦  i  +  *  (*«,)*„)  +  %(FV  +  /V+i)  (5.U) 

nv 

Weiss  and  Ascher  have  considered  the  use  of  methods  of  this  type  in  [2], [6].  The 
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collocation  methods  they  discuss  can  be  considered  as  generalizations  of  (5.11) 
and  of  the  Box  scheme,  given  by 


u 


V*-l 


-  u. 


■«JM( 


*„  +■  * 


VM 


■)(«„  +  uv¥l)  +  /  ( 


:v  + 


)• 


(5.12) 


H  A  code  based  on  those  methods  is  discussed  in  [l].  For  general  systems  (1.1) 

where  the  matrix  A(x)  is  Hermetian  or  the  equation  is  already  in  "almost”  block 
form. 


V'  = 


—An  A\z 

,  Agj  Aggj 


V 


+  F 


where 


MMxY  |  +  |Anl  +■  |A|g|  +  Mg||  +  lAgsl)  «  1 . 


centered  schemes  can  be  expected  to  behave  properly  provided  that  the  boun¬ 
dary  layers  are  properly  resolved  as  discussed  above.  However,  if  the  system 
has  not  been  blocked  beforehand,  or  if  there  are  turning  points  present  in  the 
problem,  then  in  general  we  cannot  expect  good  results  from  centered  schemes. 
The  oscillatory  nature  of  these  schemes  (see  (5.7a))  in  regions  of  the  mesh 
where  A|j4|  »  1  makes  reasonable  error  estimates  difficult  to  obtain  for  gen¬ 
eral  problems.  This  is  perhaps  best  illustrated  with  the  following  examples. 

Consider  the  system 


-1  ss  *  <  1 


(5.13) 


with  boundary  conditions 
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I 


Figure  5. 1 

V<-D  =  l.  1/(1)  =  2.  tu(-l)  =  0. 

For  e  «  1  the  solution  to  this  problem  consists  of  boundary  layers  in  the  vari- 

> 

bles  y  and  w  near  *  =  ±1  which  connect  to  the  constant  states  y  ~  0,  w  ~  0  in 
the  region  away  from  the  boundaries.  It  is  easy  to  verify,  in  fact  that 

y(x)  =  g -(»♦>)/«  +  2gl**0/«  +  o(t) 

is  the  leading  order  representation  for  y.  Since  there  is  no  non-smooth  behavior 
in  the  interior  of  the  interval,  it  might  seem  reasonable  that  this  solution  could 


i  -  lifmMii  ~  -  mini 
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be  computed  numerically  using  a  centered  scheme  and  a  nonuniform  mesh  to 
resolve  the  boundary  layers.  Figure  5.1,  which  shows  the  results  of  approximat¬ 
ing  (5.13)  with  e  =  10"a  using  the  Trapezoidal  rule,  dramatically  demonstrates 
that  this  conclusion  is  incorrect.  In  this  figure,  only  the  approximation  to  y(z) 
is  shown.  The  horizontal  lines  at  the  top  and  bottom  of  the  plot  in¬ 
dicate  the  locations  of  the  mesh  points  used  in  the  computation.  For 
scale  purposes  the  values  of  "ymin'1  and  "ymax"  in  the  legend  accompanying 
each  plot  indicate  the  locations  of  these  lines.  The  behavior  of  the  com¬ 
puted  solution  near  the  center  of  the  interval  is  clearly  unacceptable. 

One  might  suspect  that  the  behavior  near  x  =  0  is  due  to  the  potential  turn¬ 
ing  point  behavior  of  this  problem  in  that  region.  The  equations  for  y  and  v  lead 
to  the  equation 

ey"  -  *y' -  Hy  =  0  (5.U) 

for  y(x).  It  is  easy  to  verify  that  there  is  a  potential  for  nonsmooth  behavior  in 
a  neighborhood  of  size  O(VF)  near  x  =  0.  In  figure  5.2.  the  mesh  has  been 
refined  accordingly  near  z  =  0.  In  this  figure  both  y  and  tv  are  shown.  Note 
that  while  y  now  appears  smooth,  w  still  exhibits  an  unacceptable  error  near 
z  =  0.  Figure  5.3  shows  a  plot  of  the  approximation  to  w  computed  using  the 
Box  scheme  (5.12).  The  behavior  near  z  =  0  is  different  but  still  unacceptable  in 
this  case. 

It  is  possible  to  eliminate  erroneous  behavior  of  the  type  exhibited  in 
figures  5.2  and  5.3  by  using  adaptive  refinement  of  the  computational  mesh. 
However,  the  solution  will  only  become  smooth  in  the  neighborhood  of  z  =  0 
once  the  meshwidths  there  satisfy  h  =  0(c).  We  consider  this  to  be  an  unaccept- 


-36- 


able  restriction  since  in  general  turning  point  behavior  occurs  on  larger  scales 
than  0(e)  and  hence  would  require  less  refinement. 

Finally  in  figure  5.4  we  show  the  results  of  a  computation  using  the  combi¬ 
nation  of  onesided  and  centered  schemes  that  we  advocate  and  discuss  in  the 

following  sections.  One-sided  schemes  have  the  advantage  over  centered 
schemes  that  in  regions  where  h  |.A|  »  1.  they  mimic  the 'damping  behavior  of 
the  differential  equations.  This  means  that  local  errors  are  damped  out  quicUy 
by  one-sided  schemes.  In  contrast,  when  using  centered  schemes,  errors  tend 
to  be  nonlocal  due  to  the  oscillatory  behavior  of  the  methods.  As  can  be  seen 
from  the  examples  above,  this  can  result  in  significant  errors  when  solving  sys¬ 
tems  of  equations. 

Because  of  the  difficulties  of  this  type  that  can  arise  when  using  centered 
schemes,  we  have  chosen  instead  to  use  schemes  of  the  type  (5.10).  Although  it 
might  seem  that  this  involves  more  work  that  if  we  were  to  employ  centered 
schemes,  this  extra  work  is  in  fact  necessary  if  we  wish  to  guarantee  that  our 
methods  be  robust.  As  we  have  shown,  the  "cheaper”  centered  methods  can  fail 
on  problems  of  any  generality. 


-  ' 
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Figure  5. 3 


Plot  o4v»(x) 

»p$  3  0.1000E-05  Function  no.  3 

virublt  »  2  BOX  SCHEME 
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Figure  5.4 
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6.  Difference  approximations  (or  scalar  equations 

In  this  section  we  start  the  discussion  of  our  difference  approximations.  We 
divide  the  x-axis  into  subintervals  of  variable  length  hj  with  gridpoints  z,  =  0. 

v-1 

*„  =  v  ~  1*2 . N,  xn  -  c.  and  denote  by  u„  =  u(z„)  functions  defined  on 

i-o 

the  grid.  We  approximate  (2.1)  by  methods  of  the  form 

dvavuv  ¥  (1  -  cf  JavMuw+l  +  d  J v  +  (1  -  <0/wM.  (8.1) 


*.  =  y. 


We  shall  concentrate  on  two  different  methods:  the  Implicit  Euler  method  and 
the  Trapezoidal  rule  (  a  0  and  dv  «)£  respectively  ). 

We  assume  that  the  conditions  (2.2)-(2.5)  are  satisfied.  Then  it  follows  from 
the  results  of  section  2  that  the  solution  of  the  differential  equation  is  smooth 
and  we  can  obtain  error  estimates  by  standard  truncation  error  analysis.  For 
this  we  need  stability  estimates  which  we  shall  derive  now. 

Lemma  8.1:  Let  yv  i  0,  /J„  s  0  be  positive  constants  and  consider  a  grid 
function  vv  satisfying 


1  +  yv 


7*Pv 

1  +  7„' 


v  =  0.1.2 . 


or  for  0  yv  s  1 


|v^,| 


£ 


i  -  yv 

1  +  yv 


+ 


1  +7/ 


v  -  0.1.2 . 


Then 


|vj  <  max  flj  + 

0*j*  v~  l  ' 


Vs0 


with  ri  =  1/(1  +  7/)  .  ry  =  (t  -  7j )/(l  +  7j)  tor  the  first  and  second  case  respec¬ 
tively. 

Proof:  is  by  induction  on  v.  it  being  trivially  true  for  v  -  0.  For  arbitrary  u, 
both  inequalities  state  that  loy^l  s  Ty|vJ  +  (1  -  r ¥)fl¥.  with  0<tvs1.  Hence, 
using  the  induction  hypothesis, 

It^il  «  Tw(max/}w  +  Iv.IHtj)  +  (1  -  r¥)p¥ 

i <v  i<v 

=  Tymax/J,  +  (1  -  Ty)/Sy  +  !w.  I  n  Ti  ^  +  K  I  YlTt- 

j  +  v  >*v  i*  V 


Let  us  first  consider  the  case  that 


Re  a  <  -1 


(6.2) 


The  implicit  Euler  method  (dy  =  0)  can  be  written  in  the  form 


(6.3) 


(6.4) 

We  can  now  use  (6.4)  to  obtain  an  error  estimate.  Assume  thatp  &  2.  Then  y(x) 
has  two  bounded  derivatives  and 

(6.5)^^——=  +/,♦,  +  »Vm 


** 1  1  —  h¥a.¥¥ 1  1  -  A^ttvn 

ftvRea„f  |  f  ^  i 


u„ 


1  h¥n¥+ 1  1  hytiyyj  Reo 


VM 


Using  lemma  6.1  with  yv  =  |AwRea,,+,j  we  obtain  therefore 


uJ  £  max 

M  i*V 


Rea* 


i  *■  T 


7} 


where 
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k„*il  max  |y"(()l  • 

Thus  the  error  a  v  -  y v  -  uv  satisfies 


I  «J  £  max 


Ty 

Reay 


(6.8) 


This  error  estimate  is  satisfactory  if  |hyReay  |  i  1  because  in  that  case  it  follows 
that  |e„|  =  0(h*).  However  if  |hyReoty|  «  1,  then  the  method  is  not  accurate 
enough  for  our  purposes. 


We  consider  now  the  trapezoidal  rule  (d„  *  )$): 


«wi 


1  +  f^v9v*i 

1  1  -  fcMvi  * 


9v+i  -  %J  »  +  J  v+ 1). 


(6.7) 


To  estimate  |(l  +  &h„a„)/(l  -  J4h„a„«.1)|  we  have  to  distinguish  among  a  number 
of  cases. 

1)  Than  is  a  constant  as s  2  such  that  ]h„Rea„|  £  a,  )/LvReal,'.1|  £  a  for  all 
v.  In  this  case 


1  +  1  +  l 

l-J6h^avM"  1  -  \ihiAv 

where 


\ev\  =h 


b»>+i  “  au 


hiAv 

da^  dx 

1  -  V)hva-v 

1  -  | 

is  uniformly  bounded.  Let  /i„a„  =  6  +  ic ,  6.  c  real,  ]  c  |  £  p  \  b  | .  Then 

1  +  I  _  -  /—b~**~)  _  /~-^l  \~yv 

l-^a„|  V  1-6  +  J<cz+6‘  '  V  1  +  2*,  1  +  7/  K  ’ 


where 
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4 


Thus 


WK& o„|i7v=  1  +  +  68)  *  4  +  a^L  +  p*)  l^a‘<l 


l  +  Ha*“„ 

i  1 

1-7, 

1  “  HA,a,+l 

1  +  7, 

1 

'  1  +  7, 

(6.9) 


where  7„  w  yv  +  HAV6(,.  Also 


Anff 


*«■! 


/iyReavfI  9v+i 


1  —  HA,a,+i  1  ~  HA vav*\  R®®i>u 

Therefore  we  oktain  from  lemma  6. 1  for  sufficiently  small  Xthv 


I  u  J  ST  max 


H(/y  +  /i-t) 

Rea^ 

fl[l+7j 

tu.  I 


(8.10) 


where 


r-™h,trl  * (i + +p,)> + °(h> 

The  last  estimate  gives  us  again  an  error  estimate.  Assume  that  p  Jfc  3.  Then  the 
solution  of  the  differential  equation  has  three  bounded  derivatives  and  by  Taylor 
expansion  we  oktain 

"~s  H(b,+iV,m  +  “,V,)  +  H(/,m  +  /,)  +  H(+,ti  +  r„), 

* 

where 


Hl(r„*j  +  r„)|  £  rnax  |yM,(0|. 

12 


Thus 


4 


1  0 

•  B  I 


*  Tmax 

1  *j*v 


|  K(*v  +  rjJ 
j  Rea^  I' 


(6.11) 


2)  There  is  a  constant  o{  >  2  such  that  for  alt  v,  2sf/iwReav|  $a,, 
2  s  Ih^ea^il  s  a,.  This  case  can  be  reduced  to  the  previous  one  by  writing 
(6.7)  in  the  form 

mi  -  "I  1  »  V  _  gvil  1 
"  1-Avft,  l-6„  *  aw+i  1  -  c. 


where 


=  OJAvoJ'1.  c„  = 

Thus  the  same  results  as  earlier  hold,  i.e.  as  long  as  |Ji,Jlea„|  stays  bounded  we 
obtain  satisfactory  error  estimates. 

3)  |/i„Rea„|  can  become  arbitrarily  large.  In  this  case 

(1  +  jthvc l„)/(1  -  favay+x)  -  -ay/a^i 

and  the  exponential  damping  is  lost.  The  solution  of  the  homogeneous  equation 


can  grow.  We  have 


(6.11) 


vN  =  (-1 

aN 

Thus  if  ay  is  very  large  compared  with  ay,  then  tiy  is  very  large  compared  with 
vy.  This  shows  that  even  for  a  scalar  equation  the  trapezoidal  rule  need  not  be 


stable. 


The  above  error  estimates  suggest  that  we  should  use  a  combination  of  the 
trapezoidal  rule  and  the  implicit  Euler  method.  The  simplest  way  to  do  that 
would  be  to  choose  the  coefficients  dv  in  (6.1)  in  the  following  way: 

f  0  if  >  1 

if  (812) 

However  we  would  like  to  prevent  the  situation  that  we  switch  too  often  from  one 
method  to  the  other  as  a  function  of  v.  Therefore  the  are  only  allowed  to  vary 
slowly  and  we  replace  (8. 12)  by 

1)  d0  is  chosen  by  (8. 12) 

2)  For  v  >  1  we  use 

(  0  if  |/i*al(|  >  % 

If  d,_1  =  0  choose  d„  =  |j£  otherwise 

jfc  if  lhya„<  2  <813) 

//d„_lS=H  choose  d„«10  otherwise 


Assume  now  that  we  have  calculated  the  solution  of  (6. 1)  on  some  mesh 
\hv\.  Then  we  can  divide  the  interval  0  <  x  as  c  into  subintervals 
ct  <  x  ■&  c<+1,  where  ct  are  meshpoints,  i  =  0,1,  •  •  •  ,q  -1,  c„  =  0.  c,  =  c, 
such  that  on  every  subinterval  we  have  used  either  the  trapezoidal  rule  or 
the  implicit  Euler  method.  On  every  subinterval  we  can  write  down  an  error 
estimate.  These  local  error  estimates  can  be  used  to  obtain  global  error 
estimates.  Consider  an  interval  c*  <  x  :£  c1+,  and  let  y(c<)  and  u(c<)  denote 
the  solutions  at  x  -  c*  of  the  differential  equation  and  difference  approxi¬ 
mation  respectively.  Let  up  be  the  solution  of  the  difference  approximation 
with  initial  data  u ^(c4)  =  y(c<)  and  let  U#  be  the  solution  of  the  homogene¬ 
ous  difference  equation  with  u#(ct)  =  u( ct)  -up( ct).  Thus  u  =  up  +uH- 
Our  previous  results  tell  us  that 
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“i>(c<M)  =  y(c,+1)  +  tt.  it  =  0(he). 


Also 

=  V#(ct)  =  -y(ct)) 

where  |\  |  <1.  Thus  the  error  a  =  y  -  u  satisfies  the  relati»n 

«  (c<M)  =  y(citl)  -  u(c<M)  =  y  (c1+1)  -  uH(ct+1)  -  Up(q.,)  (6.14) 

=  Xt*(c<>  +  et. 

Observing  that  s(0)  =  0  we  obtain  a  linear  system  of  equatiois 


4b.  =  £ 


where 


1  0 
-A,  1 

.  0 

0  ••  •  0 

.  £  = 

e(c,) 

e(c2) 

.  £.= 

=i 

e2 

o 

•  •  •  -Vi  1 

«(ct_i) 

C9-l 

The  vector  &  represents  the  local  error.  The  global  error  is  obtained  by 
inverting  A.  In  this  particular  case  all  |X*I  <1  and  the  condition  number  is 
bounded  by  q . 

To  stress  the  interplay  between  local  and  global  error  we  remove  the 
restriction  that  Rea  £  -1.  We  allow  Rea  to  become  arbitrarily  large  and 
positive.  We  assume,  however,  that  y(r)  =  0(1)  and  slowly  varying. 
Corresponding  to  (6.12)  we  choose  the  d„  as  follows: 


0  if  |A„aJ  >  1  Rea„  <  0, 
if  !  A„a„|  s  1, 

-1  if  |A„a„|  >  1,  Rea,,  >  0. 


(6.15) 
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(In  the  same  way  as  above,  we  would  in  practice  modify  (6.17)  such  that  the 
d„  do  not  change  too  often). 

As  before  we  divide  the  interval  0  jfi  *  s  c  into'*~subintervals 
ct  £  x  <  ci  +  l  such  that  the  parameters  dv  are  constant  on  these  subinter* 
vals.  Also,  if  d„  =  )$  then  we  subdivide  ct£  x  <  ci4.l  into  subintervals 
Cy  s:  *  d  c1J+|  where  one  of  the  following  conditions  holds: 

Rea  &  -1.  |  Rea  j  <  1,  Rea  i  1. 

Without  restriction  we  can  assume  that  the  orginal  ci  are  chosen  such  that 
no  subdivision  is  necessary.  As  earlier  we  obtain  a  relation  of  the  type 
(6.14)  for  every  interval  c*  s  *  <  ct<.i  with  Rea  <  -1.  The  same  is  true  for 
intervals  with  |Rea|  <  1.  This  follows  from  well-known  error  estimates  for 
nonstiff  differential  equations.  If  Rea  >  1  then  we  use  as  initial  data 

and  integrate  the  differential  equation  in  the  direction  of  decreasing  *. 
Correspondingly  we  solve  the  difference  approximation  in  the  direction  of 
decreasing  u.  This  does  not  change  the  behavior  of  the  trapezoidal  rule  but 
the  explicit  Euler  method  d„  =  1 

(•uvM  -  uv)/h¥  =  a„u„  +  /„ 

can  be  written  as 

uv+t  fv 

—  - +•  - 

1  +  hv  av  1  +  Avai, 

and  is  the  same  as  the  implicit  Euler  method  for  decreasing  values  of  v. 
Thus  we  can  apply  our  previous  results  and  obtain 

e(c4)  =  X<M«(ciM)  +  e<,  ti  =  0(hz). 

This  shows  that  the  local  behavior  of  the  difference  approximation  is 


satisfactory  also  in  the  general  case.  However,  it  is  well  known  that  the  ini¬ 
tial  value  problem  for  (2.1)  is  not  well  posed  if  Rea  becomes  arbitrarily 
large  and  positive.  In  this  case  the  linear  system  of  equation^Jor  the  global 
error  is  not  well  conditioned.  Observe  that  the  can  be  computed  and 
therefore  also  the  condition  number  of  A  is  available. 


I 
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7.  Difference  approximations  for  diagonally  dominant  systems 

In  this  section  we  consider  systems  of  the  type  (3.2)  and  assume  that  all  the 
conditions  of  section  3  are  satisfied.  We  write  the  differential  equations  in  the 
form 

=  Ay  +  G(y )  ,  (7.1) 

where 

C(y)  =  (A  -  A)y  +  F  . 

He  approximate  (7.1)  by 

r  =  DvA{xv)uv  +  (/  —  Dv)A{x +  DVF v  +  (/  —  Dj)F „+(  (7.2) 
nv 

s  D^A(x  J)u „  +  (/  —  +  DvG(uv)  +  (/  —  Hi,)(/(U(,+i) 

Here  is  chosen  by  (8.17),  Le. 

d„(‘)  0  •  •  0 

0  0  0 

0„= . .  (7.3) 

0  . d™ 

0  if  |Avatt(x„)|  >  1,  Rea*  <  0 
dw(<)  =  H  if  \hya«{xv)\  *  1, 

1  if  |fi^aM(x„)|  >  1,  Rea«>0 


In  actual  computations  we  modify  (7.3)  in  the  same  way  as  earlier.  For  simpli¬ 
city  we  assume  that  the  dj^  =  d^  do  not  depend  on  u.  Associated  with  the  full 
system  (7.2)  is  the  diagonal  system 


By  (6.4)  and  (6.10)  we  can  find  a  constant  T  ~  2  +  p8  such  that 


|wJ<2r||(A+ A*)-'//||*  +  \vi\  +  \vff\  (7.5) 


where  we  have  used  the  notation 


||p  ||fc  =  max  |p„| . 
0*v*N 


Insiead  of  assuming  that  A  is  diagonally  dominated  we  make  a  slightly 
stronger  assumption: 

Assumption  7.1:  There  is  a  constant  6  with  0  <  6  <  1  such  that 


II (A  +  A*)-,(i4  -A)L  =fi^«/r). 


Renark:  If  we  were  only  to  use 


dW  = 


0  if  Rea*  <  0 
1  if  Rea*  >  0 


(7.6) 


then  (6.4)  would  tell  us  that  T  =  1  and  so  assumption  7.1  is  equivalent  toassum- 
ing  thatA  is  diagonally  dominated. 

In  the  same  way  as  lemma  3. 1  we  can  now  prove 

Lecuna  7.1:  If  assumption  7.1  is  valid  then  the  solutions  of  (7.2)  satisfy  the 
estimate 


l«J  <  y{,|(A+  KTl(DFv+  (I  -  D)Fv.i)\\h  +  IWI  +  lutfl).  (7.7) 


Renark:  In  the  same  way  as  for  the  continuous  problem  one  could  eAimate 
how  the  influence  of  |uj|  and  \ufj\  decreases  away  from  the  boundaries. 

Lenma  7.1  gives  us  immediately  an  error  estimate: 


-51- 


Theorem  7.1:  Assume  that  y(x)  is  a  smooth  and  bounded  solution  of  the 
differential  equation.  Then 


II V  -  till*  £  const.(A*  +  |u/  -y\ |  +  |u#  -  yff |) 


*  J 


52- 


8.  Difference  approximations  for  essentially  diagonally  dominant  systems  In 

this  section  we  consider  difference  approximations  for  the  system  (4. 1)  with 
boundary  conditions  (1.2).  They  are  of  the  form 

DvA,[zv)uv  +  (/  -  Dv)A^zv,x)uv¥X  +  DVGV  +  (/  -  DV)GV+X,  (8.1) 

fly 

-V-  --  =  +  Aa(xv.x)vv.x)  +  KEV  +  ff.M)  (8.2) 


with  boundary  conditions 


fit,, 


5. 


¥b. 


kf 


=  S 


Here  G  =  /l12v  +  F.  E  =  Azxu  +  H  and  Du  is  defined  by  (7.3).  We  assume  again 
that  D„  »  D  does  not  depend  on  v.  By  assumption  |/422h  [  «  1  and  it  is  therefore 
well-known  that  the  solution  of  (8.2)  satisfy  the  estimate 

Jv  )!*  <  RgUv,  1  +c|!£||„),  £2  ~  expdU22llo.ec). 


Lemma  7.1  tells  us  that  if  assumption  7.1  holds  then  the  solutions  of  (8.1)  satisfy 
the  estimate  (7.7)  with  F  replaced  by  G.  Using  the  assumptions  4.1  and  4.2  we 
obtain  the  analog  of  lemma  1.1: 

Lemma  8.1:  If  c  is  sufficiently  small  then  there  is  a  constant  Kg  such  that 

(IMIa  +  HujlJsAadlFlU  +  Bfflla  +  l^ol  +  l«#l  +  lv.l) 

Here  Kg  depends  only  on  Kx,  Rx  and  p. 

We  assume  again  that  the  solution  of  the  differential  equations  is  slowly 
varying  and  bounded.  Lemma  8.1  leads  then  immediately  to  the  error  estimate 

II y  -  u|U  +  !'*•■  -  v||* 

<  const.(/i2  +  \yl~ul\  +  \y([  -  uff\  +  \w9  -  vt  |).  'd'"1 


( 
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This  estimate  gives  us  an  error  estimate  for  our  problem  provided  that  the 
interval  is  sufficiently  small.  Since  again  this  cannot  be  expected  to 

be  the  case,  we  can  obtain  global  error  estimates  in  the  same  way  as  for  the 
scalar  equation  by  dividing  the  interval  up  into  subintervals.  We  divide  0  •&  x  <  c 
into  subintervals  c4  s  z  ^  ct+1,  i  =  0,1,2 . q  -1  with  the  following  properties: 

1)  D  is  constant  on  every  subinterval. 

2)  ci+1  —  ct  is  sufficiently  small  such  that  the  estimate  of  lemma  0.1  holds. 

3)  The  solution  of  the  differential  equation  is  bounded  and  slowly  varying. 

Let  y(x)  be  the  solution  of  the  differential  equation.  We  write  again 

u  =  Up  +  uH,  v  =  vP  +  Vjf  where  ( uP,vP)T  is  the  solution  of  the  difference 
approximation  with 

u/(ci)  =  v7(c<).  u}f(cui)  =  Vp(ci)=w(Ci)  (0.4) 

and  {uh,vh)t  is  the  solution  of  the  homogeneous  difference  equation  with 

uj/(c4 )  =  e'(cj,  uff( ciM)  =  e"(ctn),  vH (c4)  =  e  (c4).  (0.5) 

Here  e  =  u  —  y ,  e  =  v  -w  denotes  the  error.  Using  the  estimate  (B.3)  we 
obtain  from  (0.4) 

u#(c4)  =  y7/(Ci)  +  t(r, 

u£(c4*  i)  =  yl(cul)  +  e/.  (8.6) 

=  w(ciM)  +  c4. 

Recalling  that  the  difference  equations  are  linear  we  can  write 

tzj/(c4)  =  Lu,  u^(ctM)  =  Ll.  vff(c4  +  l)  =  L  (0.7) 

Here  L,r.  Ll ,  L  are  linear  expressions  in  e/(c4),  e77(ci+l)  and  e  (c4)  whose 
coefficients  can  be  estimated  by  A"3.  The  equations  (8.6)  and  (8.7)  give  us  n 
linear  equations 


f 


I 
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e;y(ct)  -  Ln  -  el1. 

-LJ  =  el. 

S '(c<+t)-r  =  ei 

for  every  subinterval  ct  «  x  s  ct  +  1.  There  are  g  intervals  and  n(g  +  l)  variables 
e(c4).  e  (c*).  Therefore,  using  the  boundary  conditions  we  obtain  a  linear  system 
of  n(g  +  l)  equations  in  n(g  +  l)  variables.  Again,  global  error  estimates  depend 
on  the  condition  number  of  this  system. 

Remark:  In  section  4  we  reduced  the  Solutions  of  the  differential  equations 
to  the  solution  of  a  corresponding  linear  system  of  equations.  The  two  linear 
systems  of  equations  obtained  in  that  section  and  in  this  one  need  not  be  close, 
because  we  have  not  resolved  the  potential  boundary  layers  near  x  =  Ci,cit.y  with 
an  appropriate  mesh.  Thus  in  general  the  fundamental  solutions  of  the  homo* 
geneous  differential  and  difference  equations  respectively  are  not  necessarily 
close.  However,  if  the  diagonal  elements  of  Au  satisfy  lReafl/1,,1  »  1  then  one 
can  show  using  asymptotic  expansions  that  the  corresponding  relations  are,  in 
fact,  close. 


t 
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9.  Normal  form  for  the  differential  equation 

We  will  now  discuss  how  the  general  system  (l.l)  can  be  transformed  to 
essentially  diagonally  dominant  form.  In  this  section  we  give  a  theoretical 
presentation  of  this  procedure.  The  practical  implementation  of  the  transfor¬ 
mation  differs  somewhat  from  the  discussion  in  this  section;  those  differences 
are  discussed  in  section  10. 

The  procedure  can  be  outlined  as  follows:  We  assume  that  away  from  a 
finite  number  of  turning  point  regions  the  system  is  well-behaved.  Then  the 
transformation  to  essentially  diagonally  dominant  form  is  effected  in  each  subin¬ 
terval  ofOsi  sc  through  similarity  transformations  which  put  the  matrix  A(x) 
into  an  appropriate  "blocked"  form,  and  a  stretching  of  the  independent  variable 
x  such  that  relative  to  the  basic  meshsize  h,  smoothness  requirements  similar 
to  assumptions  3. 1-3.2  and  (2.2)  are  enforced.  The  results  of  section  4  then 
guarantee  that  we  will  get  the  appropriate  error  estimates  when  the  difference 
approximation  is  applied. 

First  we  calculate  the  eigenvalues  k(x)  of  A{x)  and  divide  them  into  sets 

containing  eigenvalues  which  are  of  the  same  order  of  magnitude.  This  is 
done  in  the  following  way:  Let  K,  6  >  0  with  0  ss  Kh  «  1  be  constants.  Then 
k  G  M^0)  if  either  |*j  s  A  or  there  exists  ax  £  such  that 

|ui  -rKi|«0(ki  +  i*i).  (9-D 

By  choosing  6  sufficiently  small  we  can  guarantee  that  all  k  £  satisfy 
|xA|  «  1.  If  all  eigenvalues  belong  to  11^  then  the  construction  is  complete. 
Otherwise  let  Kt.  x2,  •  •  •  ,Kp  denote  the  eigenvalues  not  contained  in  H^0)  and  let 
I *<)  =  min  |xj.  Then  the  set  II*1*  is  formed  recursively  by  taking  k,  £  II*1*. 

Itvtp  J 

k  e  if  (Re«j)(Re«)  a  0  and  there  is  aTe  £  such  that  (9.1)  holds.  Further 
sets  are  constructed  correspondingly.  We  allow  the  number  of  sets  11^*  to 
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depend  on  z,  i.e.  as  a  function  of  x,  seta  can  split  up  and  recombine.  Therefore 
the  block-structure  can  be  a  function  of  x  as  well.  We  assume,  however,  that  we 
can  divide  the  interval  0  4  x  <  c  into  a  finite  number  of_  subintervals 
Cj  <  x  <  ctM  such  that  on  every  subinterval  the  block  structure  is  constant.  We 
will  refer  to  such  subintervals  as  " blacking  subintervals" 

The  next  step  is  to  determine  a  transformation  5  (x).  such  that 


a<*)  =  s->(xM(x)s<*)  = 


PvGO  o 

0  A-i(*) 


0 

0 


o  AGO, 


is  in  block-diagonal  form.  Here  the  eigenvalues  of  every  Aj(x)  are  exactly  the 
eigenvalues  contained  in  (counted  according  to  their  multiplicity).  We  must 
make  a  couple  of  assumptions  about  the  blocks  A  GO-  For  the  matrix  AGO  we 
require  that 


A,  l  A,  (*  >  I  «  1- 


(9-2) 


We  know  that  the  eigenvalues  of  AGO  satisfy  |«/i|  «  1.  Therefore  (9.3)  says 
that  if  we  were  to  transform  A(z)  upper  triangular  form  by  a  unitary 
transformation,  the  off-diagonal  elements  Oy  would  also  satisfy  !  OyA  |  «  1.  The 
following  shows  that  this  is  a  reasonable  assumption:  Consider  the  differential 
equation 

$-=  U'[x  )hU(z)y 


where  A  = 


0 

0 


0 


is  a  constant  matrix  and 


v  GO  = 


cosx  sinz 
l-sinx  cosz 
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is  a  unitary  matrix,  y  *  0  and  0  <  s  «  1.  The  matrix  A  is  in  upper  triangular 
form  but  the  off-diagonal  element  is  not  small.  If  we  change  variables  to 
y  =  Uy,  then  the  system  becomes  __ 


M  t 

V  = 


0  2_  +  1 
i—l  0 


f>*  ry  »v 

y  =  *y 


The  eigenvalues  of  X  are  given  by  as  =  ±  V— (y/e  +  i)  and  thus  the  solution  of  the 
equation  with  variable  coefficients  has  nothing  to  do  with  the  solution  of 


dv  /  dx  =  Aw  . 


For  the  other  blocks  we  assume  correspondingly  that 

IZj'Ajl  =  0(1).  where  E,  =  —  £  *  =  —  £  o^>.  (9.3) 

with  the  order  and  £  the  trace  of  Af.  If  these  conditions  are  not  satisfied 
then  we  choose  the  basic  meshsize  small  enough  so  that  h  \Aj  |  «  1  for  all 

Aj  that  violate  the  assumption. 

We  are  not  interested  in  highly  oscillatory  problems,  because  they  have  to 
be  treated  differently  (  see  Scheid  [5)).  Thus  we  make  the  assumption  that  the 
eigenvalues  k  of  A(x)  satisfy 

|Ime|  ip  |Re/c|  +•  C,  (9.4) 

where  p  ~  0(1)  and  0  <  C,h  «  1  are  threshold  constants.  Observe,  however, 
that  (9.4)  allows  highly  oscillatory  solutions  provided  that  the  meshsize  is 
sufficiently  small. 

We  now  describe  the  construction  of  S{x)  in  detail.  We  start  with  the  inter¬ 
val  Os  x  i  c,.  At  x  =  0  we  construct  a  unitary  transformation  U(0)  such  that 


ZT 


U*(0)A(0)U(0)  is  an  u|per  triangular  matrix  in  which  the  eigenvalues  appear  in 
the  correct  order,  i.e. 


£/'(CM(0)tf(0) 


Ar  Vr-i  1  Vo 

0  Ar-i  A- -1,0 

.  o  •••  0  A.  \ 


This  can  be  done  with  a  slightly  modified  version  of  the  usual  QR  method.  Ve 
then  determine 


2(o)  = 


i 

0 


/ 


&rJB 

‘S’r-i.O 


0  •  •  •  0  / 


such  that 


2(0)  =  5-,(0)i4(C)5(0)  = 


V 

o 


0 


0 

2r-l 


-  •  -  0 

0  0 

0  v 


S(x)  =  [/(O)^*) 


has  the  desired  block  firm.  Now  consider  the  transformed  matrix 

2(*)  =2(0)  +  B(x),  B{0)--0, 


Brr 

&r,r-l 

Bro 

B(x)  =  S~l(0)(4(*)  -A{0))S(0)  = 

Br-\,r 

Br-i,o 

Bx.t 

. 

Boo  , 

By  assumption  the  eigenvalues  of  each  block  are  well  separated  from  the  eigen¬ 
values  of  all  other  blocks.  Therefore  in  the  neighborhood  of  x  =  0  we  can  con- 
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struct  an  S(x)  such  that 


S~i(x)A(x)S(x)  = 


(*)  0 
0  0 


.  S(x)  =  S(  0)S(x) 


0  i,(x)l 


To  discuss  this  transformation  in  detail  we  neei  a  couple  ojf  lemmata. 

Lemma  9.1:  Let  An.  Agg.  E  be  pxp,  q  xg  and  pxg  matrices  respectively. 
Assume  that  the  eigenvalues  A*,!  =  1.2....,p  of  AM  are  disjoint  from  the  eigen¬ 
values  juj  j  -  1,2 . q  of  Agg.  Then  the  matrix  equation 

An  AT  -  XAgg  =  E 

has  a  unique  solution  and  there  is  a  constant  C which  depends  only  on  p,  q.  |A  |, 
|  B  |  and  mm  |  A*  -  /*,  |  such  that 

|X|  <C\E\  .  (9.5) 

Proof:  Without  restriction  we  can  assume  lhat  Atl.  Agg  are  upper  triangular; 
then  (9.5)  is  of  the  form 


Al 

®  12 

0 

A8  aZ3 

0 

*n  *12 
*a  *e2 


•  *i, 

'  *2, 


*11  *12 
*21  *22 


•  *** 


Mi 

a  i2 

*S#  ■ 

a  w 

0 

M2 

rsj 

a  g. 

0 

M, 
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For  the  first  column  of  X  we  obtain 

(Ap  —  Hi)  Xp  j  =  Bj,  j 

(*?-i  —  Ml)  2j>-l.l  =  -ap-l#2p-l  +  ep«  *— 

•  •  •  etc.  ■  ■  ■ 

Thus  the  first  column  of  X  can  be  computed  by  back-substitution  and  it  satisies 
an  estimate  of  the  type  (9.5).  The  other  columns  are  calculated  in  a  correspoid- 
ing  manner.  This  proves  the  lemma. 

Let  us  use  the  above  matrices  An  and  A&  to  form 

An  B\z 
Bg\  Agg  ' 


We  want  to  construct  a  matrix  R  such  that 


\I  R 

J  \Bzi  A 22, 

1° 

■f 


n  -  RBz\  A\\R  —  RAzz  ~  RBz\R  B \z 
Bg\  Bz\R  +  Agz 


Cu 

0 

Czi 

Cgz 

=  C 


is  a  lower  block  triangular  matrix,  i.e. 

A\\R  RAzz  RBg\R  ^  B jg  =  0. 


(97) 


Lemma  9.1  gires  us 

Lemma  t.2:  The  iteration  , 

A,,*00  -  R{n)Azz  -  R{n-l)Bz>R[n'l)  +  5, a  =  0  (98) 

converges  to  a  locally  unique  solution  of  (9.7)  provided  |£?gI|  and  Ifi^l  *re 
sufficiently  small. 
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We  shall  now  transform  C  to  block  diagonal  form. 

[«mma  9.3:  If  \Bg\\  and  \B\g\  are  sufficiently  small  then  there  is  a 
transformation  Q  such  that  v_ 


I  0 

o’ 

o 

Q 

C 1, 

o  ) 

-Q  / 

[Cn  c9z)\q  /]  =  1 

i“9^ii  +  CggQ  +  Cgj  Cgg] 

Mu  —  RBf\  0  ] 

l  0 

B2\R  +  Agg 

Proof:  If  lif^l  ,  |f?g||  are  sufficiently  small  then  the  sets  of  eigenvalues  of 
Cn  and  C22  respectively  are  still  disjoint.  Therefore  we  need  only  to  solve  the 
linear  system 


“GC11  +  CaQ  +  Cg,  =  0 

This  proves  the  lemma. 

Remark:  In  practice,  before  making  the  transformation  of  lemmata  9.1  and 
9.2,  we  diagonally  scale  the  matrix  so  that  the  off-diagonal  blocks  are  of  the 
same  order  of  magnitude:  Choosing  d  such  that  d\Bgx\  +  1  =  d'M^ul  +  1,  this 
transformation  is  given  by 


^ii 

d-‘fl,z 

(/ 

Btg 

I 

dBg\ 

An  , 

=  dI 

Bgi 

Att 

d“'/ 

This  extra  transformation  does  not  change  the  end  result  but  guarantees  that  R 
and  Q  are  of  the  same  order  of  magnitude.  The  effect  of  this  is  to  help  make 
S(z)  change  as  slowly  as  possible. 

The  results  above  can  be  used  to  construct  a  transformation  S(z)  which 
transforms  ^(x)  to  block  diagonal  form. 


* 
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Theorem  9.1:  If  14  ‘(OJtfyOOI.  t  =  1.2 . r  and  \Bqj(x)\,  j  =  0.1.2 . r  are 

sufficiently  small  then  we  can  construct  S(x)  locally  in  a  unique  way. 

Proof:  We  write  /J(x)  as  »— 


Z(x)  =  71 


pi?  a  £^22 


£  «il  «  l*r  I 

k,  «*<r> 


where 

^11  =  Ar  ^  Brr<  &  12  “  (Br.r- !•  '  '  >5ro)>  ^21  =  (.Br-i ,r*  '  ’  B^-)^  ■ 


By  assumption  |eAnl  .  I(e4i)”ll  and  i  e^zzl  are  O(l).  Also  the  eigenvalues  of 
eAn  are  well  separated  from  those  of  tAgz-  Thus  the  above  lemmata  give  us  that 
if  |e^1al  ~  |An_l5i2|  .  leia^l  ~  |A2a~li?ai|  are  sufficiently  small,  Le.  if  x  is 
sufficiently  small  then  there  is  a  unique  transformation  of  the  type 


(/  R  (/  0] 

r  +  RQ  R 

Ip  i\\q  /]  = 

Q  /, 

such  that 


Sf‘Z(z)S, 


*u  o' 

0  ^ga 

[  0  Bz\R  f  Aw 

Now,  BZIR  +■  Azz  has  the  same  properties  as  #(x)  did,  and  so  the  same  process 
can  be  applied  again  to  it.  This  proves  the  theorem.  . 

We  have  constructed  S(x)  in  a  neighborhood  of  x  =0,  but  it  is  clear  that  we 
can  continue  the  constructim  as  long  as  the  block  structure  does  not  change, 
Le.  for  Osise,.  Let  S-(«i)=  lim  S(x).  At  x  =c,  we  change  S-(c,)  to 

St(c ,)  in  the  following  way: 


t 


a 
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1)  If  two  sets  11^  merge:  S  does  not  change  (although  the  blocking  does) 

2)  If  a  set  splits  into  subsets  then  construct  a  transformation  which 
transforms  the  corresponding  block  into  block  diagonal  form.  If-necessary,  a 
permutation  matrix  can  be  applied  to  rearrange  the  blocks  according  to  the  size 
of  their  eigenvalues.  Alternatively,  5*(cj)  can  be  computed  in  the  same  way  as 
5(0). 

We  now  use  S+tej)  as  the  starting  transformation  for  the  interval 
ct£x<c 2  and  repeat  the  above  procedure  to  obtain  S(x)  at  intermediate 
points  in  that  subinterval.  In  this  way  we  determine  S(x)  everywhere  and  use  it 
to  transform  the  system  (1.1)  to 

dy  /  dx  =  A(x)y  +  H(x)y  +  G(x)  (9.9) 


with 


Ar(x)  0  •  • 

0 

0  A--1  0  ■ 

0 

0  . 

•  •  4>(*), 

H  =  -S~ldS /  dx, 
G  =  S~lF, 
y  =  S~'y 


on  every  blocking  subinterval  c4  s  x  s  c1+|. 

Ill  the  same  way  as  in  [4]  we  now  "slietch"  the  independent  vatiable.  We 
divide  each  blocking  subinterval  ct  s  x  s  c4M  into  s  i  1  "stretching  subinter- 
vals"  cy  s  x  ss  ctj+1  with  c4  =  c4#  s  c41  s  ■  ■  ■  s  cte  =  c4+I.  Suppose  we  have 
determined  cto,cu,  •  •  •  ,Cy.  Then  cij+ 4  is  determined  as  follows:  We  introduce  a 
new  variable  x  by  x  -  =  ayx  ,0s?  si  and  obtain 


dy 

dx 


=  a  A  (ax  )y  +  aH(ax)  +  aG(ax). 


(9.10) 


Here  a  =  with  0  <  a$ciH  -  c4  is  (an  approximation  to)  the  largest  value 


auch  that 
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where 


dr  I  0  •  • 

o' 

0  dr-XI  0 

0 

0  . 

•  d,I 

and  Ev  =  H(x„)uv  +  G{xv). 


The  matrix  A  is  not  in  diagonally  dominant  form.  However,  in  the  neighbor' 
hood  of  an  interior  point  x0  of  any  subinterval  we  can  find  a  constant  transfor 
mation  S  of  the  form 


S  = 


UrM  4 

Ur-xM  4-1 


'  I 


such  that  S~yA(x)S  satisfies  the  conditions  of  assumption  7.1.  Here  [/,-  is  a  uni¬ 
tary  transformation  such  that  UjAj{xa)Uj,  j  =  1,2,  •  •  •  ,r,  is  upper  triangular  for 
x  =  x0.  Lj  is  a  diagonal  scaling  such  that  LfyUjAf{x0)UjLi  satisfies  assumption 
7.1  with  }£(<5/ T)  replaced  by  I(d/r).  Then  the  smoothness  properties  of  the 
coefficients  guarantee  that  assumption  7.1  is  satisfied  in  a  whcle  neighborhood 
of  xB.  Thus  the  solution  of  the  differential  equation  satisfies  estimates  of  the 
type  (4.3)  in  the  interior  of  any  subinterval.  However,  for  interior  subintervals 
the  estimates  can  be  extended  up  to  the  boundary  provided  that 


|S*(c,)l  !5;‘(ci)l  =  0(1) 
a~y  £  5_/a+  so,  o  =  0(1)  . 


Here  a_,  5+  denote  the  stretching  factors  of  any  consecutive  subintervals.  The 
reason  the  estimates  can  be  extended  is  that  the  breakpoints  Cy  are  somewhat 
arbitrary.  We  could  move  them  a  distance  0(1/  K).  Then  the  old  breakpoints 
would  become  interior  points  in  the  new  subintervals  and  we  could  estimate  the 
derivatives.  Provided  (9.13)  holds,  these  estimates  would  not  be  destroyed  if  we 


I 
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were  to  move  the  breakpoints  back  to  the  original  position.  Finally  we  have  to 
resolve  the  boundary  layers  at  x  -  0,c .  This  is  done  as  described  at  the  end  of 
section  3. 


Jk. 


10.  Numerical  details  of  the  transformation  to  normal  form 

In  this  section  we  discuss  some  of  the  details  of  the  numerical  implementa¬ 
tion  of  the  transformation  of  a  stiff  system  to  the  normal  form  dis&ussed  in  the 
last  section.  The  blocking  and  stretching  subintervats  introduced  there  are  in 
practice  not  determined  separately  as  in  the  theoretical  discussion,  but  are 
determined  simultaneously.  For  this  reason,  we  introduce  some  new  notation. 
Starting  with  the  left  endpoint  x  =  0  =  6,  and  working  to  the  right,  we  divide  the 
interval  [O.c]  into  stretching  subintervaLs  with  endpoints  61.62.  etc.  The  block 
structure  of  the  matrix  is  monitored  as  this  is  done,  and  appropriate  points  bj 
are  designated  as  blocking  subinterval  endpoints  when  the  structure  changes. 
The  stretching  parameter  (see  (9.10))  for  the  subinterval  [6^.6yfj]  is  denoted  by 

ai  • 

An  outline  of  the  algorithm  for  the  mesh  construction  and  the  determina¬ 
tion  of  the  transformation  to  diagonally  dominant  form  follows:  We  first  deter¬ 
mine  the  eigenvalues  Kj  of  A(0).  Then  a„  is  determined  by 

-a„  minRecj  w  Kz  , 

The  stretching  parameter  an  for  the  stretching  subinterval  nearest  x  =  c  is 

determined  analogously.  Using  these  stretching  parameters  near  the  endpoints 

of  the  interval  assures  that  assumption  4.3  will  be  satisfied,  and  any  possible 

boundary  layers  will  be  resolved.  In  practice  we  construct  a  "reference  mesh” 

* 

\Xj  (jt0  with  x0  =  0,  itfn  =  c.  Xi  -  x„  =  a„  and  xjv  -  *jv-i  =  at/*  where  the  subin¬ 
tervals  [xi.Xjn]  increase  in  length  exponentially  towards  the  center  of  the  inter¬ 
val  according  to  the  rule 

x,  -  x<_,  =  min ■  (2a'0,2N~ia.lia  ) 

N*  max 

where  am#J(  is  the  maximum  length  for  a  stretching  subinterval  that  we  wish  to 


allow.  (Typically  a -  c/  10).  Then  given  a  stretching  subinterval  endpoint  by 
the  next  endpoint  by+j  is  determined  as  follows: 

A)  If  by  has  been  previously  determined  to  be  the  left  endpoinfc-of  a  blocking 
subinterval,  then  compute  S+(by)  using  QR  and  theorem  9.1  (see  remark  10.1 
below).  Let  a  =  ay_ E  =  by  +  a  be  trial  values  for  ay  and  byM  respectively. 

B)  Compute  the  eigenvalues  /c(&)  of  A(ti)  and  dete'rmine  the  sets  at 
x  =  S' .  Then 

a)  If  no  sets  M(i)  change,  go  to  C. 

b)  If  a  new  set  forms,  mark  S'  as  a  possible  blocking  subinterval  endpoint 
and  go  to  C. 

c)  If  two  sets  merge,  then  by  was  a  blocking  subinterval  endpoint.  (We  do 
not,  however,  need  to  make  a  special  computation  of  S+(6y),  since  taking 
S+(bj)  =  S-(bj)  is  acceptable  in  this  case.  The  difference  in  the  treatment 
of  this  subinterval  is  that  S(S)  will  be  computed  by  updating  S^by)  taking 
the  new  block  structure  into  account.)  Go  to  C. 

C)  Compute  the  transformation  matrix  S (S')  by  updating  5 (by)  using 
theorem  9.1  (see  remarks  10.1  and  10.2).  We  need  this  updated  version  of  S 
even  if  x  =  21  is  a  blocking  subinterval  endpoint 

D)  Now  compute  the  left-hand  sides  of  the  tests  in  (9.11).  (The  actual  imple¬ 

mentation  of  these  tests  is  discussed  in  remark  10.3).  During  the  determination 
of  byM  from  by  it  is  possible  to  return  to  this  step  (D)  several  times.  Using  the 
reference  mesh,  we  first  determine  5  ■  where  by  e  [xt,xi+1].  The 

action  taken  can  be  different  in  each  case: 

Dl)  (First  time):  If  any  test  fails,  replace  3  with  3/V2  (i.e.  decrease  the 
length  of  the  stretching  subinterval)  and  try  again  (go  to  B).  If  no  tests  fait, 
then  replace  3  with  min(vT?  a,o)  (  increase  3)  and  go  to  B. 

D2)  (Second  time):  If  no  tests  failed  under  Dl  but  any  test  fails  this  time, 
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set  ay  =  a/  VS  and  go  to  E  (below).  If  no  test  failed  this  time,  replace  a 
with  min(V2a,5)  and  go  to  B).  If  any  tests  failed  under  D1  and  any  test  fails 
this  time,  replace  a  bya/  V2  and  go  to  R 

D3)  (Third  time):  If  no  tests  fail,  set  ay  =  a  and  go  to  E.  If  any  test  fails  at 
this  point,  but  no  test  failed  under  D2).  set  ay  =  a/ V2  and  go  to  E.  If  any 
test  fails  at  this  point  and  any  test  failed  under  Dl.  t(ien  special  action  must 
be  taken,  because  if  we  were  to  decrease  a  any  further,  we  would  have 
a.j/ aj- 1  ^  }f.  We  would  like  to  avoid  this  situation  for  the  reasons  discussed 
at  the  end  of  the  last  section  (with  a  =  2).  First,  however,  we  must  deter¬ 
mine  how  small  a  actually  needs  to  be  near  bj.  So  replace  a  with  a/ 2  and 
go  to  B. 

D4)  (Fourth  and  succeeding  times):  If  any  test  fails,  replace  a  with  a/  2  and 
goto  B.  If  no  tests  fail,  then  we  have  a  value  for  a  which  represents  the 
proper  stretching  near  x  -  bf.  In  order  to  use  this,  however,  we  have  to 
redistribute  the  previous  endpoints  so  that  there  is  a  smooth 

exponential  grading  of  subintervals  into  the  region  around  *  =  6^.  We  do 
this  by  first  looking  for  the  minimum  value  of  i  with 

2‘S  *  2a*., 

where  k  is  determined  by  first  calculating 

x=bJ-‘S^2l  . 

and  then  finding  k  such  that  x  E  Then  the  stretching  subinterval 

endpoints  bk+i,bk+z,  ■  ,bj  are  replaced  by 

fej+a/Z*-1  .  I  =  1.2.  ,i 

=  \bh  +  a/2‘-‘  ,  I=i  +  1,  ,m 


where  m  is  chosen  such  that  <  6^  £  6t+(n.  Steps  E  and  F  must  now 

be  redone  for  all  of  these  corrected  subintervals  [bt.bt+i]. 
l  =  k  +  1,  •  •  •  ,k  +m.  Then  set  j  =  k  +m  and  go  to  A. 

E)  The  subinterval  endpoint  has  now  been  determined,  i.e.  6j  +  1  =  &. 
Stretch  the  interval  +  to  0  fi  5  s  1,  Put  down  a  uniform  mesh 
with  meshwidth  K  where  K  is  a  meshwidth  that  would  be  oonsidered  appropriate 
for  the  resolution  of  a  smooth  function  on  the  interval  0  i  x  <  1.  (This  is  some¬ 
what  vague;  typically  K  ~  1/  10).  The  meshpoints  in  the  original  variable  x  are 
given  by 

zW  =  bj  +  i/fiaj  v  =  0,1 1  /K. 

FO  Now  compute  the  transformation  matrices  S(x^l),  v  -  1,2 . l/£  by 

updating  S+(bj)  using  theorem  9.1  (again  see  remarks  10.1  and  10.2).  The 
difference  approximation  can  now  be  written  down  for  the  mesh  intervals 

].  1/  =  0,1,2 . 1/Ji-l.  Suppressing  the  superscript  | j],  the  difference 

approximation  is  given  by 

V  —  V  v  —  DV(hVAy+\  Sv+  1  (5„H  —  S  V))v  „+j 

+  (/  -  Dv)(hvAv  -  S;'  (5^1  -  Sv))v  „  (10.1) 

+  hvDvSv+iF +  h „(/  -  Py)S v  1 F y 

where  Sv  =  S(xt^)  and  y  v  js  an  approximation  to  y(zv).  In  practice  we  have 
found  it  more  convenient  to  make  our  computations  in  terms  of  the  original 

x 

(untransformed)  variables  y(x„)  which  we  approximate  with  v„.  The  difference 
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approximation  then  becomes 

Sy+iVmi  —  Sv*  vv  —  +  Sv+\  —  5,, 1  )v  „+ , 

+  (/  -  Dv)(hvS~lAv  +  5JV‘i  -  Sk1  )vv  (10.2) 

+  hvDvSv+\Fv+\  +  Aw(/  ~  Dv)SvlFv 

where  A*  =  A(xJ^).  This  is  done  because  it  is  usually  the  original  variables  that 
we  are  interested  in.  One  should  note,  however,  that  if  lSj  +  { S’J'1  |  is  not  of 
reasonable  size,  one  can  expect  the  system  (10.1)  to  be  somewhat  better  condi¬ 
tioned  that  (10.2).  The  transformed  variables  y  are  in  some  sense  the  "correct" 
variables  for  the  problem  since  they  have  been  scaled  in  such  away  that  we  can 
obtain  estimates  for  the  system. 

At  this  point  we  can  now  increment  j  and  return  to  step  A. 

Remarks  10.1:  When  computing  the  transformation  S*( by)  at  the  left  end¬ 
point  of  a  blocking  subinterval,  we  first  transform  A{b})  to  upper  triangular  form 
using  the  QR  method.  If  the  eigenvalues  do  not  appear  in  the  correct  order  on 
the  diagonal  of  the  transformed  matrix,  in  practice  we  repeat  the  QR  iteration 
using  the  (now  known)  eigenvalues  as  shifts  in  the  order  in  which  we  wish  them 
to  appear. 

The  resulting  matrix  is  then  transformed  to  block  diagonal  form  using  lem¬ 
mata  9.2  and  9.3.  Note  that  if  the  eigenvalue  sets  11^  are  well  separated  then 
the  iteration  (9.8)  can  be  replaced  with 

+  R<n-l)BtlR{n-"  -  Bti,  (10.3) 

i.e.  we  only  need  to  invert  An.  (This  remark  also  applies  when  updating  S  later 
on  ). 

Note  also  that  the  ofi-diagonal  blocks  need  not  be  completely  eliminated  as 
any  0(1)  blocks  can  be  absorbed  into  the  matrix  H  of  (9.9). 


-72- 


Remark  10.8  As  long  as  a  blocking  subinterval  endpoint  does  not  lie 
between  two  points  Xj,Xj+1,  S (2j«.i)  can  always  be  computed  from  S(xj)  by 
updating  S(zj)  using  the  blocking  technique  of  lemmata  9.2,9.3.  In  practice  the 
iteration  (9.8)  is  replaced  by 

+  S12  =  o  (10.4) 


where 


*n=  U[AuVy,  Xgg=U^2gUi 

are  upper  triangular  (  and  Ug  are  unitary  and  are  determined  by  QR).  Here 
3\z=  U[B\zUz  and  Sgi-  UgBziU\.  R  is  then  computed  from  3  using 
R  =  UjftUg.  This  simplifies  the  computation  since  in  particular  the  iteration 
(10.4)  only  involves  the  solution  of  systems  of  the  form  (9.6).  The  only  exception 
is  if  the  eigenvalue  sets  11^  are  well  enough  separated  so  that  (10.3)  can  be  used 
in  place  of  (9.8). 

Remark  10.3:  In  practice  we  only  compute  the  tests  (9.11)  with  p  =  1. 
Although  this  means  that  we  do  not  necessarily  get  the  smoothness  required  for 
good  error  estimates,  our  experience  has  been  that  we  always  obtain  satisfac¬ 
tory  results.  Also  taking  p  =  1  means  that,  the  difference  approximations  to 
(9.11)  will  only  involve  two  adjacent  subinterval  endpoints  bj,  which  simplifies 
the  algorithm  considerably.  The  difference  approximations  for  the  smoothness 
tests  are  described  as  follows:  Suppose  we  want  to  assure  that  the  function  /(x) 
satisfies 


(max(|a/(ox  )|,l))_larf/(ax  )/ dx  «  K. 
In  practice  we  replace  this  with  the  test 

5|/(6j  +  a)  -/(6,)| 

i  +  ai/(6,)i 


(10.5) 


(10.8) 
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for  some  appropriate  value  of  ft »  K,  where  bj  is  the  previously  determined 
subinterval  endpoint  and  a  is  the  stretching  parameter  that  we  are  testing. 
(10.8)  is  obtained  by  stretching  [bj.bj  +  a]  to  [0,1]  and  then  ..replacing  the 
derivative  of  (10.5)  with  a  divided  difference  over  the  whole  interval  0  s  z  s  1. 
The  denominator  of  (10.5)  has  been  replaced  by  a  sum  because  it  is  cheaper  to 
compute  than  the  maximum  but  gives  approximately  the  same  effect. 

We  have  computed  several  examples  to  test  our  procedures,  and  we  present 
three  of  them  here.  Each  example  was  chosen  to  test  a  different  aspect  of  the 
transformation  to  normal  form.  For  the  first  example  we  consider  the  system 


€ 

fc  +  3*2 


-lszsl 


with  y(-l)  =  1,  |/(1)  =  2.  This  problem  has  turning  point  behavior  near 
x  =  -V2.0.V2,  and  thus  we  are  testing  the  aspect  of  the  mesh  construction  that 
looks  at  the  smoothness  of  the  coefficients  in  order  to  determine  the  proper 
stretching.  Figure  (10.1)  shows  the  results  of  a  computation  with  e  =  10-8.  Only 
the  approximation  to  y(z)  is  shown.  As  with  the  earlier  plots,  the  horizontal 
lines  above  and  below  the  plots  are  used  to  indicate  the  locations  of  the  the 
meshpoints. 

The  second  example  (figure  10.2)  shows  an  example  with  both  boundary 
layers  and  a  possible  turning  point  at  *  =0.  The  system  is  given  by 


d_y 

dx  v 


X _ J_ 

e  e 

V 

0 

y 

y(-l)=  1.  V (1)  =  2.  e  =  10-3. 


Note  that  the  region  of  the  possible  turning  point  has  been  refined  even  though 
the  solution  is  smooth  there.  In  order  that  our  code  be  robust  we  have  chosen 
to  "resolve"  all  possible  unsmooth  behavior  even  though  in  cases  such  as  this 


one  it  might  be  possible  to  tell  a  priori  thst  the  solution  will  be  smooth  there. 


In  the  third  example  (figures  10.3,  10.4),  we  test  the  feature  of  the  mesh 
refinement  algorithm  that  resolves  possible  highly  oscillatory  behavior.  The  sys¬ 
tem  we  consider  is  given  by 


»(-D  *  l.y(l)  =2.  e  =  10-V 


In  this  example,  the  eigenvalues  of  the  natrix  become  complex  in  a  neighbor- 

1. 

hood  of  x  =  0  of  width  0(e*).  The  mesh  in  this  region  is  thus  determined  by 
imposing  the  fourth  of  conditions  (9. 11)  Another  interesting  feature  of  this 
example  is  that  it  is  not  particularly  well-posed.  The  solution  becomes  very 
large  near  the  left  boundary  (||y(x)||-  =®.9xl0#).  However,  the  method  is  able 
to  handle  the  situation  quite  well.  Figure  10.4  is  a  magnification  of  the  region 
near  x  =  0.  This  shows  the  oscillations  that  occur  near  x  =  0,  and  demonstrates 
that  the  mesh  has  been  properly  scaled  toresolve  these  oscillations. 
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E  dd(y)/dxdx  ♦  <x*3  -  x/2)  dy/dx  -  y  »  8 
E*  epsilon  «  0.1000E-04  No.  ol  Mshpoint*  *  124 
ynin  «  -0.8543E-01  ymx  *  0 .1871E+01 


Figure  10.1 
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£  dd<y)/dxdx  ♦  <x‘2)  dy/dx  ♦  y  «  0 
E  *  tptilofl  *  0.1000E-03  No.  o<  anhpoinU  *  229 
y»in  *  -Q.8270E+0?  yiux  ■  0.B912EMB 


Figure  10.  S 


E  dd<y)/dxdx  «  <x‘2)  dy/dx  ty«0 

E  3  epsilon  3  O.iOOOE-03  No.  of  iMshpoints  3  130 
Min  3  -0.1212E*00  mix  3  0.4471E«00 

Min  3  -0.9179E+06  Mix  3  0.2473E+07 


Figure  10.4 
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11.  Solution  of  nonlinear  systems 

It  is  relatively  straight-forward  to  apply  our  methods  to  nonlinear  problems 
of  the  form 

■^V(x)  =  f{y(z))  ■  (11.1) 

where  for  simplicity  we  specify  linear  boundary  conditions  (12).  Here  y  and  / 
are  vector  functions  of  dimension  n.  As  in  [3],  we  solve  (11  1), (1.2)  using  a  func¬ 
tional  Newton  iteration  technique:  We  linearize  (11.1)  about  a  previous  guess  or 
approximate  solution  yn(x),  obtaining 

=  ^fiyn{x))y(x)  +  /(yn(x))  -  ~yn(z)  .  (11.2) 

Here  y(x)  =  y"*1  -yn  is  the  correction  to  the  guess  yn.  and  |^(yn(x))  =  A(x) 

is  the  Jacobian  matrix  of  / .  Letting  F  =  /  (yB(x))  -  -^-yn(x)  we  see  that  (1 1.2) 

is  in  the  same  form  as  (1.1)  and  so  we  can  apply  the  methods  discussed  earlier 
in  this  paper  for  linear  problems.  We  emphasize  here  that  the  linearization  is 
done  before  the  method  is  applied.  This  is  in  contrast  to  the  usual  approach  for 
solving  noninear  ODEs,  which  is  to  first  apply  the  difference  method,  and  then 
use  a  Newton  iteration  to  solve  the  resulting  system  of  nonlinear  algebraic  equa¬ 
tions.  The  reason  for  this  is  that  we  can  only  expect  our  method  to  apply  to 
linear  differential  equations,  and  so  to  be  certain  that  it  works,  we  must  first 
reduce  the  nonlinear  ODEs  to  linear  ones. 

We  have  used  this  approach  to  solve  a  simple  test  problem,  given  by 


ty"  +  jiy2)' -y  =0,  -lsis  l,  y(-i)  =  l.  y(i)  =  2. 
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As  in  [3]  we  replace  this  with  the  2x2  system  of  equations 

ty'  =  -jyz 
v1  =  y 


(11.3) 


and  then  linearize  this  system.  The  matrix  transformations  and  mesh  construc¬ 
tion  are  the  same  as  before.  The  difference  approximation  is  essentially  the 
same  as  (10.2)  except  that  the  term  dyn/dx  is  grouped  with  the  term  dyAix 
when  making  the  difference  approximation.  In  terms  of  vv  and  wjf ,  the  approxi¬ 
mations  to  y(xv )  and  yn(xv)  respectively,  and  letting  Fv  =  f  {yn{xv)),  we  add  the 
terms 


-  Sv-* )  -  S&  )y^x  +  ((/  -  DV)(S^\  -  S;1 )  +  S'1 )y? 

to  the  right-hand  side  of  (10.2)  to  obtain  our  difference  approximation. 

Figure  11.1  shows  the  result  of  a  numerical  computation  of  the  solution  of 
(11.3)  for  e  =  10-3.  We  started  with  an  initial  guess  for  y  of  a  straight  line 
between  the  boundary  values  and  continued  in  e  with  the  values  e  =  .1,  .03,  .015, 
.0075,  .004.  .001.  Convergence  was  obtained  to  a  tolerance  of  less  than  10“3  in 
the  maximum  norm  at  each  step 


J. 
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